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Abstract: We study chemical reactions with complex mechanisms under two assumptions: 
(i) intermediates are present in small amounts (this is the quasi-steady-state hypothesis or 
QSS) and (ii) they are in equilibrium relations with substrates (this is the quasiequilibrium 
hypothesis or QE). Under these assumptions, we prove the generalized mass action law 
together with the basic relations between kinetic factors, which are sufficient for the 
positivity of the entropy production but hold even without microreversibility, when the 
detailed balance is not applicable. Even though QE and QSS produce useful approximations 
by themselves, only the combination of these assumptions can render the possibility beyond 
the "rarefied gas" limit or the "molecular chaos" hypotheses. We do not use any a priori form 
of the kinetic law for the chemical reactions and describe their equilibria by thermodynamic 
relations. The transformations of the intermediate compounds can be described by the 
Markov kinetics because of their low density (low density of elementary events). This 
combination of assumptions was introduced by Michaelis and Menten in 1913. In 1952, 
Stueckelberg used the same assumptions for the gas kinetics and produced the remarkable 
semi-detailed balance relations between collision rates in the Boltzmann equation that 
are weaker than the detailed balance conditions but are still sufficient for the Boltzmann 
//-theorem to be valid. Our results are obtained within the Michaelis-Menten-Stueckelbeg 
conceptual framework. 

Keywords: chemical kinetics; Lyapunov function; entropy; quasiequilibrium; detailed 
balance; complex balance 
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1. Introduction 

1.1. Main Asymptotic Ideas in Chemical Kinetics 

There are several essentially different approaches to asymptotic and scale separation in kinetics, and 
each of them has its own area of applicability. 

In chemical kinetics various fundamental ideas about asymptotical analysis were developed [1]: 
Quasieqiulibrium asymptotic (QE), quasi steady-state asymptotic (QSS), lumping, and the idea of 
limiting step. 

Most of the works on nonequilibrium thermodynamics deal with the QE approximations and 
corrections to them, or with applications of these approximations (with or without corrections). There 
are two basic formulation of the QE approximation: The thermodynamic approach, based on entropy 
maximum, or the kinetic formulation, based on selection of fast reversible reactions. The very first 
use of the entropy maximization dates back to the classical work of Gibbs [2], but it was first 
claimed for a principle of informational statistical thermodynamics by Jaynes [3]. A very general 
discussion of the maximum entropy principle with applications to dissipative kinetics is given in the 
review [4]. Corrections of QE approximation with applications to physical and chemical kinetics were 
developed [5,6]. 

QSS was proposed by Bodenstein in 1913 [7], and the important Michaelis and Menten work [8] 
was published simultaneously. It appears that no kinetic theory of catalysis is possible without QSS. 
This method was elaborated into an important tool for the analysis of chemical reaction mechanism 
and kinetics [9-11]. The classical QSS is based on the relative smallness of concentrations of some 
of the "active" reagents (radicals, substrate-enzyme complexes or active components on the catalyst 
surface) [12-14]. 

Lumping analysis aims to combine reagents into "quasicomponents" for dimension 
reduction [15,17-19]. Wei and Prater [16] demonstrated that for (pseudo)monomolecular 
systems there exist linear combinations of concentrations which evolve in 
time independently. These linear combinations (quasicomponents) correspond 

to the left eigenvectors of the kinetic matrix: If IK = \l then 
d(/, c)/dt = (I, c)A, where the standard inner product (I, c) is the concentration of a quasicomponent. 
They also demonstrated how to find these quasicomponents in a properly organized experiment. 

This observation gave rise to a question: How to lump components into proper quasicomponents 
to guarantee the autonomous dynamics of the quasicomponents with appropriate accuracy? Wei 
and Kuo studied conditions for exact [15] and approximate [17] lumping in monomolecular and 
pseudomonomolecular systems. They demonstrated that under certain conditions a large monomolecular 
system could be well-modelled by a lower-order system. 

More recently, sensitivity analysis and Lie group approach were applied to lumping analysis [18,19], 
and more general nonlinear forms of lumped concentrations were used (for example, concentration of 
quasicomponents could be a rational function of c). 

Lumping analysis was placed in the linear systems theory and the relationships between lumpability 
and the concepts of observability, controllability and minimal realization were demonstrated [20]. 
The lumping procedures were considered also as efficient techniques leading to nonstiff systems and 
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demonstrated the efficiency of the developed algorithm on kinetic models of atmospheric chemistry 
[21]. An optimal lumping problem can be formulated in the framework of a mixed integer nonlinear 
programming (MINLP) and can be efficiently solved with a stochastic optimization method [22] . 

The concept of limiting step gives the limit simplification: The whole network behaves as a single 
step. This is the most popular approach for model simplification in chemical kinetics and in many areas 
beyond kinetics. In the form of a bottleneck approach this approximation is very popular from traffic 
management to computer programming and communication networks. Recently, the concept of the 
limiting step has been extended to the asymptotology of multiscale reaction networks [23,24]. 

In this paper, we focus on the combination of the QE approximation with the QSS approach. 

1.2. The Structure of the Paper 

Almost thirty years ago one of us published a book [25] with Chapter 3 entitled "Quasiequilibrium 
and Entropy Maximum". A research program was formulated there, and now we are in the position to 
analyze the achievements of these three decades and formulate the main results, both theoretical and 
applied, and the unsolved problems. In this paper, we start this work and combine a presentation of 
theory and application of the QE approximation in physical and chemical kinetics with exposition of 
some new results. 

We start from the formal description of the general idea of QE and its possible extensions. In 
Section 2, we briefly introduce main notations and some general formulas for exclusion of fast variables 
by the QE approximation. 

In Section 3, we present the history of the QE and the classical confusion between the QE 
and the quasi steady state (QSS) approximation. Another surprising confusion is that the famous 
Michaelis-Menten kinetics was not proposed by Michaelis and Menten in 1913 [8] but by Briggs 
and Haldane [12] in 1925. It is more important that Michaelis and Menten proposed another 
approximation that is very useful in general theoretical constructions. We described this approximation 
for general kinetic systems. Roughly speaking, this approximation states that any reaction goes through 
transformation of fast intermediate complexes (compounds), which (i) are in equilibrium with the input 
reagents and (ii) exist in a very small amount. 

One of the most important benefits from this approach is the exclusion of nonlinear kinetic laws and 
reaction rate constants for nonlinear reactions. The nonlinear reactions transform into the reactions of the 
compounds production. They are in a fast equilibrium and the equilibrium is ruled by thermodynamics. 
For example, when Michaelis and Menten discuss the production of the enzyme-substrate complex ES 
from enzyme E and substrate S, they do not discuss reaction rates. These rates may be unknown. They 
just assume that the reaction E + S ^ ES is in equilibrium. Briggs and Haldane involved this reaction 
into the kinetic model. Their approach is more general than the Michaelis-Menten approximation but 
for the Briggs and Haldane model we need more information, not only the equilibrium of the reaction 
E + S ^ ES but also its rates and constants. 

When compounds undergo transformations in a linear first order kinetics, there is no need to include 
interactions between them because they are present in very small amounts in the same volume, and their 
concentrations are also small. (By the way, this argument is not applicable to the heterogeneous catalytic 
reactions. Although the intermediates are in both small amounts and in a small volume, i.e., in the 
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surface layer, the concentration of the intermediates is not small, and their interaction does not vanish 
when their amount decreases [33]. Therefore, kinetics of intermediates in heterogeneous catalysis may 
be nonlinear and demonstrate bifurcations, oscillations and other complex behavior.) 

In 1952, Stueckelberg [26] used similar approach in his seminal paper "if -theorem and unitarity of 
the S-matrix". He studied elastic collisions of particles as the quasi-chemical reactions 

v + w — > v' + w' 

(v, w, v', w' are velocities of particles) and demonstrated that for the Boltzmann equation the linear 
Markov kinetics of the intermediate compounds results in the special relations for the kinetic coefficients. 
These relations are sufficient for the //-theorem, which was originally proved by Boltzmann under the 
stronger assumption of reversibility of collisions [27]. 

First, the idea of such relations was proposed by Boltzmann as an answer to the Lorentz objections 
against Boltzmann's proof of the //-theorem. Lorentz stated the nonexistence of inverse collisions for 
polyatomic molecules. Boltzmann did not object to this argument but proposed the "cyclic balance" 
condition, which means balancing in cycles of transitions between states Si — > S 2 — > ■ ■ ■ — > S n — > 
Si. Almost 100 years later, Cercignani and Lampis [28] demonstrated that the Lorenz arguments 
are wrong and the new Boltzmann relations are not needed for the polyatomic molecules under the 
microreversibility conditions. The detailed balance conditions should hold. 

Nevertheless, Boltzmann's idea is very seminal. It was studied further by Heitler [29] and Coester [30] 
and the results are sometimes cited as the "Heitler-Coestler theorem of semi-detailed balance". In 1952, 
Stueckelberg [26] proved these conditions for the Boltzmann equation. For the micro-description he 
used the S'-matrix representation, which is in this case equivalent for the Markov microkinetics (see 
also [31]). 

Later, these relations for the chemical mass action kinetics were rediscovered and called the complex 
balance conditions [51,63]. We generalize the Michaelis-Menten-Stueckelberg approach and study in 
Section 5 the general kinetics with fast intermediates present in small amount. In Subsection 5.7 the big 
Michaelis-Menten-Stueckelberg theorem is formulated as the overall result of the previous analysis. 

Before this general theory, we introduce the formalism of the QE approximation with all the necessary 
notations and examples for chemical kinetics in Section 4. 

The result of the general kinetics of systems with intermediate compounds can be used wider than 
this specific model of an elementary reaction: The intermediate complexes with fast equilibria and the 
Markov kinetics can be considered as the "construction staging" for general kinetics. In Section 6, we 
delete the construction staging and start from the general forms of the obtained kinetic equations as from 
the basic laws. We study the relations between the general kinetic law and the thermodynamic condition 
of the positivity of the entropy production. 

Sometimes the kinetics equations may not respect thermodynamics from the beginning. To repair this 
discrepancy, deformation of the entropy may help. In Section 7, we show when is it possible to deform 
the entropy by adding a linear function to provide agreement between given kinetic equations and the 
deformed thermodynamics. As a particular case, we got the "deficiency zero theorem". 

The classical formulation of the principle of detailed balance deals not with the thermodynamic and 
global forms we use but just with equilibria: In equilibrium each process must be equilibrated with 
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its reverse process. In Section 7, we demonstrate also that for the general kinetic law the existence of a 
point of detailed balance is equivalent to the existence of such a linear deformation of the entropy that the 
global detailed balance conditions (Equation (87) below) hold. Analogously, the existence of a point of 
complex balance is equivalent to the global condition of complex balance after some linear deformation 
of the entropy. 

1.3. Main Results: One Asymptotic and Two Theorems 

Let us follow the ideas of Michaelis-Menten and Stueckelberg and introduce the asymptotic theory of 
reaction rates. Let the list of the components Aj be given. The mechanism of reaction is the list of the 
elementary reactions represented by their stoichiometric equations: 



The linear combinations J2i a pi^-i an d ^ /3 P i A are the complexes. For each complex ^\ y^Ai from 
the reaction mechanism we introduce an intermediate auxiliary state, a compound B y Each elementary 
reaction is represented in the form of the "2n-tail scheme" (Figure 1) with two intermediate compounds: 



There are two main assumptions in the Michaelis-Menten-Stueckelberg asymptotic: 

• The compounds are in fast equilibrium with the corresponding input reagents (QE); 

• They exist in very small concentrations compared to other components (QSS). 

The smallness of the concentration of the compounds implies that they (i) have the perfect 
thermodynamic functions (entropy, internal energy and free energy) and (ii) the rates of the reactions 
B.[ — > Bj are linear functions of their concentrations. 

One of the most important benefits from this approach is the exclusion of the nonlinear reaction 
kinetics: They are in fast equilibrium and equilibrium is ruled by thermodynamics. 

Under the given smallness assumptions, the reaction rates r p for the elementary reactions have a 
special form of the generalized mass action law (see Equation (74) below): 




(1) 




(2) 



Figure 1. A 2n-tail scheme of an extended elementary reaction. 




r p = y2 p exp(a p ,/2) 
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where ip p > is the kinetic factor and exp(a p , /2) is the Boltzmann factor. Here and further in the text 
(a p , fi) = a pifa is the standard inner product, exp( , ) is the exponential of the value of the inner 
product and fa are chemical potentials fi divided on RT. 

For the prefect chemical systems, fa = ln(cj/c*), where q is the concentration of Ai and c* > are 
the positive equilibrium concentrations. For different values of the conservation laws there are different 
positive equilibria. The positive equilibrium c* is one of them and it is not important which one is it. At 
this point, fa = 0, hence, the kinetic factor for the perfect systems is just the equilibrium value of the 
rate of the elementary reaction at the equilibrium point c*\ cp p = r p (c*). 

The linear kinetics of the compound reactions £>j — > Bj implies the remarkable identity for the 
reaction rates, the complex balance condition (Equation (89) below) 

^2 exp(/2, a p ) = ^2^P P exp(/2, f3 p ) 
p p 

for all admissible values of fi and given ip which may vary independently. For other and more convenient 
forms of this condition see Equation (91) in Section 6. The complex balance condition is sufficient 
for the positivity of the entropy production (for decrease of the free energy under isothermal isochoric 
conditions). The general formula for the reaction rate together with the complex balance conditions and 
the positivity of the entropy production form the Michaelis-Menten-Stueckelberg theorem (Section 5.7). 
The detailed balance conditions (Equation (87) below), 

for all p, are more restrictive than the complex balance conditions. For the perfect systems, the detailed 
balance condition takes the standard form: r p (c*) = r~(c*). 

We study also some other, less restrictive sufficient conditions for accordance between 
thermodynamics and kinetics. For example, we demonstrate that the G-inequality (Equation (92) below) 

^<p p exp(/2,a p ) > ^<p p exp(/2,/3 p ) 
p p 

is sufficient for the entropy growth and, at the same time, weaker than the condition of complex balance. 

If the reaction rates have the form of the generalized mass action law but do not satisfy the sufficient 
condition of the positivity of the entropy production, the situation may be improved by the deformation of 
the entropy via addition of a linear function. Such a deformation is always possible for the zero deficiency 
systems. Let q be the number of different complexes in the reaction mechanism, d be the number of the 
connected components in the digraph of the transitions between compounds (vertices are compounds and 
edges are reactions). To exclude some degenerated cases a hypothesis of weak reversibility is accepted: 
For any two vertices Bi and Bj, the existence of an oriented path from B { to Bj implies the existence of 
an oriented path from Bj to B^ 

Deficiency of the system is [63] 

q — d — rankr > 

where T = (7^) = {(3^ — cc^) is the stoichiometric matrix. If the system has zero deficiency then 
the entropy production becomes positive after the deformation of the entropy via addition of a linear 
function. The deficiency zero theorem in this form is proved in Section 7.3. 
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Interrelations between the Michaelis-Menten-Stueckelberg asymptotic and the transition state theory 
(which is also referred to as the "activated-complex theory", "absolute-rate theory", and "theory of 
absolute reaction rates") are very intriguing. This theory was developed in 1935 by Eyring [35] and 
by Evans and Polanyi [36]. 

Basic ideas behind the transition state theory are [37]: 

• The activated complexes are in a quasi-equilibrium with the reactant molecules; 

• Rates of the reactions are studied by studying the activated complexes at the saddle point of a 
potential energy surface. 

The similarity is obvious but in the Michaelis-Menten-Stueckelberg asymptotic an elementary 
reaction is represented by a couple of compounds with the Markov kinetics of transitions between them 
versus one transition state, which moves along the "reaction coordinate", in the transition state theory. 
This is not exactly the same approach (for example, the theory of absolute reaction rates uses the detailed 
balance conditions and does not produce anything similar to the complex balance). 

Important technical tools for the analysis of the Michaelis-Menten-Stueckelberg asymptotic are the 
theorem about preservation of the entropy production in the QE approximation (Section 2 and Appendix 
1) and the Morimoto //-theorem for the Markov chains (Appendix 2). 

2. QE and Preservation of Entropy Production 

In this section we introduce informally the QE approximation and the important theorem about the 
preservation of entropy production in this approximation. In Appendix 1, this approximation and the 
theorem are presented with more formal details. 

Let us consider a system in a domain U of a real vector space E given by differential equations 

The QE approximation for (3) uses two basic entities: Entropy and slow variables. 

Entropy S is an increasing concave Lyapunov function for (3) with non-degenerated Hessian 

d 2 S/dxidxf 

dS 

In this approach, the increase of entropy in time is exploited (the Second Law in the form (4)). 

The slow variables M are defined as some differentiable functions of variables x: M = m{x). Here 
we assume that these functions are linear. More general nonlinear theory was developed in [38,39] with 
applications to the Boltzmann equation and polymer physics. Selection of the slow variables implies a 
hypothesis about separation of fast and slow motion. The slow variables (almost) do not change during 
the fast motion. After some initial time, the fast variables with high accuracy are functions of the slow 
variables: We can write x m x* M . 

The QE approximation defines the functions x* M as solutions to the following MaxEnt 
optimization problem: 

S(x) — > max subject to m(x) = M (5) 
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The reasoning behind this approximation is simple: During the fast initial layer motion, entropy increases 
and M almost does not change. Therefore, it is natural to assume that x* M is close to the solution to the 
MaxEnt optimization problem (5). Further x* u denotes a solution to the MaxEnt problem. 

A solution to (5), x* M , is the QE state, the set of the QE states x* M , parameterized by the values of the 
slow variables M is the QE manifold, the corresponding value of entropy 

S*{M) = S(x* M ) (6) 

is the QE entropy and the equation for the slow variables 

m(F(x* M )) (7) 



dt 

represents the QE dynamics. 

The crucial property of the QE dynamics is the preservation of entropy production. 

Theorem about preservation of entropy production. Let us calculate dS*(M)/dt at point M 
according to the QE dynamics (7) and find dS(x)/dt at point x = x* M due to the initial system (3). 
The results always coincide: 

dS*(M) dS(x 



dt dt 



(8) 



The left hand side in (8) is computed due to the QE approximation (7) and the right hand side 
corresponds to the initial system (3). The sketch of the proof is given in Appendix 1. 

The preservation of the entropy production leads to the preservation of the type of dynamics: If for 
the initial system (3) entropy production is non-negative, dS/dt > 0, then for the QE approximation (7) 
the production of the QE entropy is also non-negative, dS* /dt > 0. 

In addition, if for the initial system (dS/dt)\ x = if and only if F(x) = then the same property 
holds in the QE approximation. 

3. The Classics and the Classical Confusion 

3.1. The Asymptotic of Fast Reactions 

It is difficult to find who introduced the QE approximation. It was impossible before the works of 
Boltzmann and Gibbs, and it became very well known after the works of Jaynes [3]. 

Chemical kinetics has been a source for model reduction ideas for decades. The ideas of QE appear 
there very naturally: Fast reactions go to their equilibrium and, after that, remain almost equilibrium all 
the time. The general formalization of this idea looks as follows. The kinetic equation has the form 

^L=K sl {N) + -K fs {N) (9) 
dt e 

Here N is the vector of composition with components Ni > 0, K s i corresponds to the slow reactions, 
Kf s corresponds to fast reaction and e > is a small number. The system of fast reactions has the linear 
conservation laws b t (N) = £\ b^Nf. k(K fs (N)) = 0. 
The fast subsystem 
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tends to a stable positive equilibrium iV* for any positive initial state N(0) and this equilibrium is a 
function of the values of the linear conservation laws fy(iV(0)). In the plane bi(N) = bi(N(0)) the 
equilibrium is asymptotically stable and exponentially attractive. 

Vector b(N) = (bi(N)) is the vector of slow variables and the QE approximation is 

f t = b(K sl (N*(b)) (10) 

In chemical kinetics, equilibria can be described by conditional entropy maximum (or conditional 
extremum of other thermodynamic potentials). Therefore, for these cases we can apply the formalism of 
the quasiequilibrium approximation. The thermodynamic Lyapunov functions serve as tools for stability 
analysis and for model reduction [40]. 

The QE approximation, the asymptotic of fast reactions, is well known in chemical kinetics. Another 
very important approximation was invented in chemical kinetics as well. It is the Quasi Steady State 
(QSS) approximation. QSS was proposed in [7] and was elaborated into an important tool for analysis 
of chemical reaction mechanisms and kinetics [9-11]. The classical QSS is based on the relative 
smallness of concentrations of some of "active" reagents (radicals, substrate-enzyme complexes or active 
components on the catalyst surface) [13,14]. In the enzyme kinetics, its invention was traditionally 
connected to the so-called Michaelis-Menten kinetics. 

3.2. QSS and the Briggs-Haldane Asymptotic 

Perhaps the first very clear explanation of the QSS was given by Briggs and Haldane in 1925 [12]. 
Briggs and Haldane consider the simplest enzyme reaction S + E ±=> SE — >■ P + E and mention that 
the total concentration of enzyme ([E] + [SE]) is "negligibly small" compared with the concentration 
of substrate [S]. After that they conclude that -^[SE] is "negligible" compared with -^[S] and ^[P] and 
produce the now famous 'Michaelis-Menten' formula, which was unknown to Michaelis and Menten: 
h[E}[S] = + k 2 )[ES] or 

where the "Michaelis-Menten constant" is 



K 



M 



There is plenty of misleading comments in later publications about QSS. Two most important 
confusions are: 

• Enzymes (or catalysts, or radicals) participate in fast reactions and, hence, relax faster than 
substrates or stable components. This is obviously wrong for many QSS systems: For example, 
in the Michaelis-Menten system all reactions include enzyme together with substrate or product. 
There are no separate fast reactions for enzyme without substrate or product. 

• Concentrations of intermediates are constant because in QSS we equate their time derivatives 
to zero. In general case, this is also wrong: We equate the kinetic expressions for some time 
derivatives to zero, indeed, but this just exploits the fact that the time derivatives of intermediates 
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concentrations are small together with their values, but not obligatory zero. If we accept QSS then 
these derivatives are not zero as well: To prove this we can just differentiate the Michaelis-Menten 
formula (11) and find that [ES] in QSS is almost constant when [S] 3> Km, this is an additional 
condition, different from the Briggs-Haldane condition [E] + [AE] <C [S] (for more details 
see [1,14,33] and a simple detailed case study [41]). 

After a simple transformation of variables the QSS smallness of concentration transforms into a 
separation of time scales in a standard singular perturbation form (see, for example [33,34]). Let us 
demonstrate this on the traditional Michaelis-Menten system: 

M = - kl [S}[E}+k„ 1 [SE] 

^ = k 1 [S][E]-(k^ 1 + k 2 )[SE] W 
[E\ + [SE\ = e = const, [S] + [P] = s = const 

This is a homogeneous system with the isochoric (fixed volume) conditions for which we write the 
equations. The Briggs-Haldane condition is e C s. Let us use dimensionless variables x = [S]/s, 
£ = [SE]/c: 

= ~ sk Ml - + fc_i£ 
JL d< (13) 

■± = 8k lX (l-Z)-(k- 1 + k 2 )Z 

To obtain the standard singularly perturbed system with the small parameter at the derivative, we need 
to change the time scale. This means that when e — > the reaction goes proportionally slower and to 
study this limit properly we have to adjust the time scale: dr = -dt: 



dx 

— = -shx{i - + fc-i£ 

s dr 



(14) 



For small e/s, the second equation is a fast subsystem. According to this fast equation, for a given 
constant x, the variable £ relaxes to 

sx 

£qss 



K M + sx 

exponentially, as exp(— (sk\x + k_ x + ^2)^)- Therefore, the classical singular perturbation theory 
based on the Tikhonov theorem [42,43] can be applied to the system in the form (14) and the QSS 
approximation is applicable even on an infinite time interval [44]. This transformation of variables and 
introduction of slow time is a standard procedure for rigorous proof of QSS validity in catalysis [33], 
enzyme kinetics [45] and other areas of kinetics and chemical engineering [13]. 

It is worth to mention that the smallness of parameter e/s can be easily controlled in experiments, 
whereas the time derivatives, transformation rates and many other quantities just appear as a result of 
kinetics and cannot be controlled directly. 
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3. 3. The Michaelis and Menten Asymptotic 

QSS is not QE but the classical work of Michaelis and Menten [8] was done on the intersection of 
QSS and QE. After the brilliantly clear work of Briggs and Haldane, the name "Michaelis -Menten" 
was attached to the Briggs and Haldane equation and the original work of Michaelis and Menten was 
considered as an important particular case of this approach, an approximation with additional and not 
necessary assumptions of QE. From our point of view, the Michaelis -Menten work includes more and 
may give rise to an important general class of kinetic models. 

Michaelis and Menten studied the "fermentative splitting of cane sugar". They introduced 
three "compounds": The sucrose-ferment combination, the fructose-ferment combination and the 
glucose-ferment combination. The fundamental assumption of their work was "that the rate of 
breakdown at any moment is proportional to the concentration of the sucrose-invertase compound". 

They started from the assumption that at any moment according to the mass action law 

[S l ][E]=K l [S l E] (15) 

where [Si] is the concentration of the ith sugar (here, i = for sucrose, 1 for fructose and 2 for glucose), 
[E] is the concentration of the free invertase and A'j is the ith equilibrium constant. 

For simplification, they use the assumption that the concentration of any sugar in question in free state 
is practically equal to that of the total sugar in question. 

Finally, they obtain 

K (l + q[P]) + [Sq] 
where e = [E] + [P] = [Si] = [S 2 ] and q = £ + 

Of course, this formula may be considered as a particular case of the Briggs -Haldane formula (11) 
if we take ^> k 2 in (11) (i.e., the equilibration S + E <=? SE is much faster than the reaction 
SE —7- P + E) and assume that q = in (16) (i.e., fructose-ferment combination and glucose-ferment 
combination are practically absent). 

This is the truth but may be not the complete truth. The Michaelis-Menten approach with many 
compounds which are present in small amounts and satisfy the QE assumption (15) is a seed of the 
general kinetic theory for perfect and non-perfect mixtures. 

4. Chemical Kinetics and QE Approximation 

4.1. Stoichiometric Algebra and Kinetic Equations 

In this section, we introduce the basic notations of the chemical kinetics formalism. For more details 
see, for example, [33]. 

The list of components is a finite set of symbols A±, . . . , A n . 

A reaction mechanism is a finite set of the stoichiometric equations of elementary reactions: 

^2 a p*Ai -> ^2 PpiAi (17) 

i i 

where p = l,...,m is the reaction number and the stoichiometric coefficients a P i,[3 P i are 
nonnegative integers. 
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A stoichiometric vector "f p of the reaction (17) is a n-dimensional vector with coordinates 

l P i = (3 P i - a pi (18) 

that is, "gain minus loss" in the pth elementary reaction. 

A nonnegative extensive variable iVj, the amount of Ai, corresponds to each component. We call the 
vector N with coordinates Ni "the composition vector". The concentration of Ai is an intensive variable 
q = Ni/V, where V > is the volume. The vector c = N/V with coordinates q is the vector of 
concentrations. 

A non-negative intensive quantity, r p , the reaction rate, corresponds to each reaction (17). The kinetic 
equations in the absence of external fluxes are 

p 

If the volume is not constant then equations for concentrations include V and have different form (this is 
typical for the combustion reactions, for example). 

For perfect systems and not so fast reactions, the reaction rates are functions of concentrations and 
temperature given by the mass action law for the dependance on concentrations and by the generalized 
Arrhenius equation for the dependance on temperature T. 

The mass action law states: 

r p {c,T) = k p {T)\{c a r (20) 

i 

where k p (T) is the reaction rate constant. 
The generalized Arrhenius equation is: 

k p {T)=A p exp(^e X p(-^ (21) 

where R = 8.314472 K ^ ol is the universal, or ideal gas constant, E ap is the activation energy, S ap is 
the activation entropy (i.e., E ap — TS ap is the activation free energy), A p is the constant pre-exponential 
factor. Some authors neglect the S ap term because it may be less important than the activation energy, 
but it is necessary to stress that without this term it may be impossible to reconcile the kinetic equations 
with the classical thermodynamics. 

In general, the constants for different reactions are not independent. They are connected by various 
conditions that follow from thermodynamics (the second law, the entropy growth for isolated systems) or 
microreversibility assumption (the detailed balance and the Onsager reciprocal relations). In Section 6.2 
we discuss these conditions in more general settings. 

For nonideal systems, more general kinetic law is needed. In Section 5 we produce such a general 
law following the ideas of the original Michaelis and Menten paper (this is not the same as the famous 
"Michaelis-Menten kinetics"). For this work we need a general formalism of QE approximation for 
chemical kinetics. 
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4. 2. Formalism of QE Approximation for Chemical Kinetics 
4.2.1. QE Manifold 

In this section, we describe the general formalism of the QE for chemical kinetics following [34]. 

The general construction of the quasi-equilibrium manifold gives the following procedure. First, let us 
consider the chemical reactions in a constant volume under the isothermal conditions. The free energy 
F(N,T) = Vf(c,T) should decrease due to reactions. In the space of concentrations, one defines a 
subspace of fast motions L. It should be spanned by the stoichiometric vectors of fast reactions. 

Slow coordinates are linear functions that annulate L. These functions form a subspace in the space 
of linear functions on the concentration space. Dimension of this space is s — n — dim L. It is necessary 
to choose any basis in this subspace. We can use for this purpose a basis bj in L- 1 , an orthogonal 
complement to L and define the basic functionals as bj(N) = (bj, N). 

The description of the QE manifold is very simple in the Legendre transform. The chemical potentials 
are partial derivatives 

dF(N,T) df(c,T) 
^ = ^— = ^c— (22) 
Let us use fa as new coordinates. In these new coordinates (the "conjugated coordinates"), the QE 
manifold is just an orthogonal complement to L. This subspace, L L , is defined by equations 

^ Hili = for any 7 G L (23) 

i 

It is sufficient to take in (23) not all 7 6 L but only elements from a basis in L. In this case, we get the 
system of n — dim L linear equations of the form (23) and their solution does not cause any difficulty. 
For the actual computations, one requires the inversion from \i to c. 

It is worth to mention that the problems of the selection of the slow variables and of the description 
of the QE manifold in the conjugated variables can be considered as the same problem of description of 
the orthogonal complement, L L . 

To finalize the construction of the QE approximation, we should find for any given values of slow 
variables (and of conservation laws) bi the corresponding point on the QE manifold. This means that we 
have to solve the system of equations for c: 

b(N)=b; (/i(c,T), 7p ) = (24) 

where b is the vector of slow variables, [i is the vector of chemical potentials and vectors 7 P form a basis 
in L. After that, we have the QE dependence cqe(6) and for any admissible value of b we can find all 
the reaction rates and calculate b. 

Unfortunately, the system (24) can be solved analytically only in some special cases. In general case, 
we have to solve it numerically. For this purpose, it may be convenient to keep the optimization statement 
of the problem: F — > min subject to given b. There exists plenty of methods of convex optimization for 
solution of this problem. 

The standard toy example gives us a fast dissociation reaction. Let a homogeneous reaction 
mechanism include a fast reaction of the form A + B =f± AB. We can easily find the QE approximation 
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for this fast reaction. The slow variables are the quantities b\ = N A — N B and b 2 = N A + N B + N c 
which do not change in this reaction. Let the chemical potentials be ha/RT = lnc^ + hao, Hb/RT = 
lncB+/i_B , ^ab/RT = In cab + At abo- This corresponds to the free energy F = VRT J2 i q (In Q+/i i0 ), 
the correspondent free entropy (the Massieu-Planck potential) is —F/T. The stoichiometric vector is 
7 = ( — 1, — 1, 1) and the equations (24) take the form 

h b 2 cab „ 

ca~ c b = — , c A + c B + cab = tj , = K 

V V c A c B 

where K is the equilibrium constant K = exp(/XAo + f^Bo — £*abo)- 

From these equations we get the expressions for the QE concentrations: 



CA{b\M) 



l&i l Ifih i\ 2 h + b 2 




IV K V V 2 V K KV 



cb(6i, b 2 ) = c A (bi, b 2 ) - y , cab(&i, ^2) = &1 ^ &2 - 2ca(&i, 62) 
The QE free entropy is the value of the free entropy at this point, c(&i, 62)- 

4.2.2. QE in Traditional MM System 

Let us return to the simplest homogeneous enzyme reaction E + S ^ ES P + S, the traditional 
Michaelis-Menten System (12) (it is simpler than the system studied by Michaelis and Menten [8]). 
Let us assume that the reaction E + S ^ ES is fast. This means that both k\ and k_i include large 
parameters: k\ = -k±, k_i = For small e, we will apply the QE approximation. Only three 

components participate in the fast reaction, A\ = S, A 2 = E, A 3 = ES. For analysis of the QE 
manifold we do not need to involve other components. 

The stoichiometric vector of the fast reaction is 7 = (—1,-1,1). The space L is one-dimensional 
and its basis is this vector 7. The space L L is two-dimensional and one of the convenient bases is 
bi = (1, 0, 1), b 2 = (0, 1, 1). The corresponding slow variables are &i(iV) = JVi + N 3 , b 2 {N) = N 2 + JV 3 . 
The first slow variable is the sum of the free substrate and the substrate captured in the enzyme- substrate 
complex. The second of them is the conserved quantity, the total amount of enzyme. 

The equation for the QE manifold is (15): fc 1 c 1 c 2 = /c_ic 3 or %£| = p- because fcic^c^ = fc-iCg, 
where c* = c*(T) > are the so-called standard equilibrium values and for perfect systems 
tH = RT\n{ Ci /c*), F = RTV ^ c^cjc*) - 1). 

Let us fix the slow variables and find ci j2 ,3. Equations (24) turn into 

ci + c 3 = 61 , c 2 + c 3 = b 2 , k x c x c 2 = fc_ic 3 

Here we change dynamic variables from N to c because this is a homogeneous system with 
constant volume. 

If we use ci = &i — c 3 and c 2 = b 2 — c 3 then we obtain a quadratic equation for c 3 : 



fciCg - (kibi + k x b 2 + k_i)c 3 + fci&i& 2 = 



(25) 
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Therefore, 



Abib 2 



The sign "— " is selected to provide positivity of all Cj. This choice provides also the proper asymptotic: 
c 3 — > if any of foj — > 0. For other c lj2 we should use c\ — b\ — c 3 and c 2 = b 2 — c 3 . 
The time derivatives of concentrations are: 

Ci = -A4C1C2 + fc_]C 3 + WinC 1 ] 11 - 17 out Ci 

C 2 = -A4C1C2 + (fc-1 + fc 2 )c 3 + t'i n C 2 n - f out C 2 

(26) 

c 3 = A;icic 2 - (A;_i + k 2 )c 3 + v in c 3 n - v out c 3 
c 4 = &2C 3 + Vi n C™ - v out c 4 

here we added external flux with input and output velocities (per unite volume) v- m and v out and input 
concentrations c m . This is done to stress that the QE approximation holds also for a system with fluxes 
if the fast equilibrium subsystem is fast enough. The input and output velocities are the same for all 
components because the system is homogeneous. 
The slow system is 

b\ = ci + c 3 = -k 2 c 3 + v m H? - v out bi 



v in b 2 r 



b 2 = c 2 + c 3 

Q = k 2 C 3 + WinC™ - 

where &i n = cf + c 3 n , 6 2 n = c 2 n + c 3 n . 

Now, we should use the expression for c 3 (»i, 6 2 ): 

1 



- U ou t&2 
^out c 4 



(27) 



61 



c 4 =fc 2 



61 + 62 



fc_i 




&1 + &2 



^out02 




4&1&2 



(28) 



111^4 



It is obvious here that in the reduced system (28) there exists one reaction from the lumped component 
with concentration bi (the total amount of substrate in free state and in the substrate-enzyme complex) 
into the component (product) with concentration C4. The rate of this reaction is k 2 c(bib 2 ). The 
lumped component with concentration b 2 (the total amount of the enzyme in free state and in the 
substrate-enzyme complex) affects the reaction rate but does not change in the reaction. 

Let us use for simplification of this system the assumption of the substrate excess (we follow the logic 
of the original Michaelis and Menten paper [8]): 



[S] > [SE] 



i.e., 



bi > c 3 



Under this assumption, the quadratic equation (25) transforms into 

k-i 



b 2 

1 + r 
t»i 



c 3 



■ C3 



(29) 



(30) 
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and in this approximation 

6i + h + — 

(compare to (16) and (11): This equation includes an additional term b 2 in denominator because we did 
not assume formally anything about the smallness of b 2 in (29)). 
After this simplification, the QE slow equations (28) take the form 

• h 2 b 2 b x . , in , 

bi = t— + VjJ) x - w out 6i 

bi+b 2 + k -^ 

b 2 = v in b™ - v out b 2 (32) 
h + b 2 + k -^ 

This is the typical form in the reduced equations for catalytic reactions: Nominator in the reaction rate 
corresponds to the "brutto reaction" S + E ->> P + E [33,49]. 

4.2.3. Heterogeneous Catalytic Reaction 

For the second example, let us assume equilibrium with respect to the adsorption in the CO on Pt 
oxidation: 

CO+Pt^PtCO; 2 +2Pt^2PtO 

(for detailed discussion of the modeling of CO on Pt oxidation, this "Mona Liza" of catalysis, we address 
readers to [33]). The list of components involved in these 2 reactions is: A\ = CO, A 2 = 2 , A 3 = Pt, 
A A = PtO, A 5 = PtCO (C0 2 does not participate in adsorption and may be excluded at this point). 

Subspace L is two-dimensional. It is spanned by the stoichiometric vectors, 71 = (—1, 0, —1, 0, 1), 
72 = (0,-1,-2,2,0). 

The orthogonal complement to L is a three-dimensional subspace spanned by vectors (0, 2, 0, 1, 0), 
(1, 0, 0, 0, 1), (0, 0, 1, 1, 1). This basis is not orthonormal but convenient because of integer coordinates. 
The corresponding slow variables are 

61 = 2N 2 + N 4 = 2N Q2 + N PtQ 

b 2 = N, + iV 5 = N co + iVptco (33) 

b 3 = N 3 + iV 4 + N 5 = N Pt + N Pt0 + iVp tC o 

For heterogeneous systems, caution is needed in transition between N and c variables because there are 
two "volumes" and we cannot put in (33) q instead of iV^: iV gas = V gas c gas but iV sm .f = V^rfCsurf, where 
where V gas is the volume of gas, V sur f is the area of surface. 

There is a law of conservation of the catalyst: iVpt + iVp t o + ^Vptco = ^3 = const. Therefore, we 
have two non-trivial dynamical slow variables, 61 and b 2 . They have a very clear sense: b\ is the amount 
of atoms of oxygen accumulated in 2 and PtO and b 2 is the amount of atoms of carbon accumulated in 
CO and PtCO. 

The free energy for the perfect heterogeneous system has the form 

F = V gas RT c i ( ln R) " l) + VsurfRT c * ( ln ( ~* ) " 1 ) ( 34 ) 

^ lg as V \ C i/ / A . surf \ \ C iJ J 



Entropy 2011, 1 3 



982 



where q are the corresponding concentrations and c* = c*(T) > are the so-called standard equilibrium 
values. (The QE free energy is the value of the free energy at the QE point.) 
From the expression (34) we get the chemical potentials of the perfect mixture 

li^RTlnl^) (35) 



c ■ 

The QE manifold in the conjugated variables is given by equations: 

-fJ>i - /U 3 + /i 5 = ; -ju 2 - 2fi 3 + 2/i 4 = 
It is trivial to resolve these equations with respect to /i 3 4, for example: 



or with the standard equilibria: 



1 

A*4 = + A*3 ; ^5 = ^1 + ^3 



C4 c 3 /C 2 C 5 Ci c 3 



c 4 c 3 V C 2 c 5 c l c 3 



or in the kinetic form (we assume that the kinetic constants are in accordance with thermodynamics and 
all these forms are equivalent): 

kiCic 3 = /c_ic 5 , k 2 c 2 c\ = k_ 2 c 2 4 (36) 
The next task is to solve the system of equations: 

fcicic 3 = k-!C 5 , k 2 c 2 cl = k_ 2 c\ , 2V g3S c 2 + V snri c A = b x , 

VgasCl + V; ur fC5 = b 2 , V; ur f(c 3 + C 4 + C 5 ) =63 

This is a system of five equations with respect to five unknown variables, 01,2,3,4,5. We have to solve them 
and use the solution for calculation of reaction rates in the QE equations for the slow variables. Let us 
construct these equations first, and then return to (37). 

We assume the adsorption (the Langmuir-Hinshelwood) mechanism of CO oxidation (the numbers in 
parentheses are used below for the numeration of the reaction rate constants): 

(±1) CO+Pt^PtCO 

(±2) 2 +2Pt^2PtO (38) 
(3) PtO+PtCO^C0 2 +2Pt 

The kinetic equations for this system (including the flux in the gas phase) is 



CO 


Ni 


= V am{ (-kicic 3 + fc_ic 5 ) 


+ V gas (v in cf - UoutCi) 


o 2 


N 2 


= V suvi (-k 2 c 2 c 2 3 + k_ 2 c\) 


+ ^gasl^inC™ — w out c 2) 


Pt 


N 3 


= Kurf(-A:iCiC 3 + /c_iC 5 • 


- 2fc 2 c 2 Cg + 2A;_ 2 C4 + 2k 3 c A c 5 ) 


PtO 


Na 


= V smi (2k 2 c 2 cl - 2k_ 2 c\ 


— k 3 c^c<~) 


PtCO 


N 5 


= ^surf(^lClC 3 - /c_iC 5 - 




co 2 


N 6 


— ^s Ur f/c 3 C4C5 + V^as (Wi n Cg n — V ut c 6) 



(39) 
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Here v- m and v out are the flux rates (per unit volume). 
For the slow variables this equation gives: 

hi = 2N 2 + N4 = -V SUT{ k 3 c 4 c 5 + 2^ gas (i; in C2 n - v out c 2 ) 

b 2 = Ni + N 5 = -V snri k 3 C4;C 5 + VgasOinC 1 / 1 - U out Ci) 

63 = N 3 + iV 4 + N 5 = 

N 6 = V BVx{ k 3 c 4 c b + V gas (v- m c^ - v out c 6 ) 

This system looks quite simple. Only one reaction, 

PtO+PtCO^C0 2 +2Pt (41) 

is visible. If we know expressions for c 3 5(6) then this reaction rate is also known. In addition, to work 
with the rates of fluxes, the expressions for c lj2 (b) are needed. 

The system of equations (37) is explicitly solvable but the result is quite cumbersome. Therefore, let 
us consider its simplification without explicit analytic solution. We assume the following smallness: 

61 > N 4 , b 2 > N 5 (42) 

Together with this smallness assumptions equations (37) give: 

h 



c 3 



v i 1 -i- fc i b 2 i A Js b i 

Vsurf \ J- T fc l y gas -T W 2 fc _ 2 y gas 



1 fc 2 61 63 
C4 - A/ o~TZv~~ 71 n Z z \ (43) 



c 5 



2 fc_ 2 \/ gas y A + + Jlp-jL-) 

fci 6 2 63 



In this approximation, we have for the reaction (41) rate 



. , h 1 k 2 Vhb 2 b\ 

r = k 3 c A c b = k 3 - 



1 k 2 61 

2 fc_2 Vo-as 



This expression gives the closure for the slow QE equations (40). 

4.2.3. Discussion of the QE procedure for Chemical Kinetics 

We finalize here the illustration of the general QE procedure for chemical kinetics. As we can see, the 
simple analytic description of the QE approximation is available when the fast reactions have no joint 
reagents. In general case, we need either a numerical solver for (24) or some additional hypotheses about 
smallness. Michaelis and Menten used, in addition to the QE approach, the hypothesis about smallness of 
the amount of intermediate complexes. This is the typical QSS hypothesis. The QE approximation was 
modified and further developed by many authors. In particular, a computational optimization approach 
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for the numerical approximation of slow attracting manifolds based on entropy-related and geometric 
extremum principles for reaction trajectories was developed [47]. 

Of course, validity of all the simplification hypotheses is a crucial question. For example, for the 
CO oxidation, if we accept the hypothesis about the quasiequilibrium adsorption then we get a simple 
dynamics which monotonically tends to the steady state. The state of the surface is unambiguously 
presented as a continuous function of the gas composition. The pure QSS hypothesis results for the 
Langmuir-Hinshelwood reaction mechanism (38) without quasiequilibrium adsorption in bifurcations 
and the multiplicity of steady states [33]. The problem of validity of simplifications cannot be 
solved as a purely theoretical question without the knowledge of kinetic constants or some additional 
experimental data. 

5. General Kinetics with Fast Intermediates Present in Small Amount 

5. 1 . Stoichiometry of Complexes 

In this Section, we return to the very general reaction network. 

Let us call all the formal sums that participate in the stoichiometric equations (17), the complexes. The 
set of complexes for a given reaction mechanism (17) is 0i, . . . , Q q . The number of complexes q < 2m 
(two complexes per elementary reaction, as the maximum) and it is possible that q < 2m because some 
complexes may coincide for different reactions. 

A complex O; is a formal sum Oj = YTj=i u ij^-j = i u h where ^ is a vector with coordinates z/jy. 

Each elementary reaction (17) may be represented in the form 0~ — > 0+ where 0^ are the 
complexes which correspond to the right and the left sides (17). The whole mechanism is naturally 
represented as a digraph of transformation of complexes: Vertices are complexes and edges are reactions. 
This graph gives a convenient tool for the reaction representation and is often called the "reaction graph". 

Let us consider a simple example: 18 elementary reactions (9 pairs of mutually reverse reactions) 
from the hydrogen combustion mechanism (see, for example, [48]). 



H + 2 ^ O + OH; O + H 2 ^ H + OH; 

OH + H 2 ^H + H 2 0; O + H 2 ^ 20H; 

H0 2 + H^H 2 + 2 ; H0 2 + H^20H; 



(44) 



H + OH + M ^ H 2 + M; H + 2 + M ^ H0 2 + M; 
H 2 2 + H ^ H 2 + H0 2 



There are 16 different complexes here: 







i 



H + 2 , 2 = O + OH, 3 = O + H 2 , 4 = H + OH, 
OH + H 2 , 6 = H + H 2 0, 7 = O + H 2 0, 8 = 20H, 



©5 







9 



H0 2 + H, 0io = H 2 + 2 , On = H + OH + M, 







12 



H 2 + M, 



13 — 



H + 2 + M, 0i4 = H0 2 + M, 







15 



H 2 2 + H, 0i6 = H 2 + H0 2 
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The reaction set (44) can be represented as 

Oi ^ e 2 , e 3 ^ e 4 , e 5 ^ e 6 , e 7 - e 8 ^ © 9 ^ e 10 , 

011 ^ 012, ©13 ^ ©14, ©15 ^ ©16 

We can see that this digraph of transformation of complexes has a very simple structure: There are five 
isolated pairs of complexes and one connected group of four complexes. 

5. 2. Stoichiometry of Compounds 

For each complex Qj we introduce an additional component Bj, an intermediate compound and B^ 
are those compounds Bj (1 < j < q), which correspond to the right and left sides of reaction (17). 

We call these components "compounds" following the English translation of the original 
Michaelis-Menten paper [8] and keep "complexes" for the formal linear combinations Bj. 

An extended reaction mechanism includes two types of reactions: Equilibration between a complex 
and its compound (q reactions, one for each complex) 

&j ^ Bj (45) 

and transformation of compounds B~ — > (m reactions, one for each elementary reaction from (17). 
So, instead of the reaction (17) we can write 

i i 

Of course, if the input or output complexes coincide for two reactions then the corresponding 
equilibration reactions also coincide. 

It is useful to visualize the reaction scheme. In Figure 1 we represent the 2n-tail scheme of an 
elementary reaction sequence (46) which is an extension of the elementary reaction (17). 

The reactions between compounds may have several channels (Figure 2): One complex may transform 
to several other complexes. 

The reaction mechanism is a set of multichannel transformations (Figure 2) for all input complexes. 
In Figure 2 we grouped together the reactions with the same input complex. Another representation of 
the reaction mechanism is based on the grouping of reactions with the same output complex. Below, in 
the description of the complex balance condition, we use both representations. 

The extended list of components includes n + q components: n initial species Ai and q compounds 
Bj. The corresponding composition vector N® is a direct sum of two vectors, the composition vector 
for initial species, N, with coordinates jVj (i — 1, . . . , n) and the composition vector for compounds, T, 
with coordinates Tj (j = 1, . . . , q): N® = N ® T. 

The space of composition vectors E is a direct sum of n-dimensional Ea and g-dimensional E B : 
E = E A ® E B . 

For concentrations of Ai we use the notation q and for concentrations of Bj we use <^-. 
The stoichiometric vectors for reactions Qj ^ Bj (45) are direct sums: g J = — Vj © ej, where 6j is 
the jth standard basis vector of the space R q = E B , the coordinates of ej are eji = Sjf. 

g 1 = {-Vjl, ~Vj2, • • • , -Vj n , 0, ■ ■ ■ , 0, 1, 0, . . . , 0) (47) 

i 
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The stoichiometric vectors of equilibration reactions (45) are linearly independent because there exists 
exactly one vector for each I. 

The stoichiometric vectors 7 j7 of reactions Bj — > Bi belong entirely to Eb- They have jth coordinate 
— 1, /th coordinate +1 and other coordinates are zeros. 

To exclude some degenerated cases a hypothesis of weak reversibility is accepted. Let us consider 
a digraph with vertices 0j and edges, which correspond to reactions from (17). The system is weakly 
reversible if for any two vertices 9; and Qj, the existence of an oriented path from Oj to Qj implies the 
existence of an oriented path from Bj to 6j. 

Of course, this weak reversibility property is equivalent to weak reversibility of the reaction network 
between compounds Bj. 

Figure 2. A multichannel view on the complex transformation. The hidden reactions 
between compounds are included in an oval S. 



5.3. Energy, Entropy and Equilibria of Compounds 

In this section, we define the free energy of the system. The basic hypothesis is that the compounds 
are the small admixtures to the system, that is, the amount of compounds Bj is much smaller than 
amount of initial components A { . Following this hypothesis, we neglect the energy of interaction between 
compounds, which is quadratic in their concentrations because in the low density limit we can neglect 
the correlations between particles if the potential of their interactions decay sufficiently fast when the 
distance between particles goes to oo [50]. We take the energy of their interaction with A4 in the linear 
approximation, and use the perfect entropy for Bi. These standard assumptions for a small admixtures 
give for the free energy: 





(48) 
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Let us introduce the standard equilibrium concentrations for Bj. Due to the Boltzmann distribution 
(exp(-u/BT)) and formula (48) 

<-(c,r) = iexpf-^l (49) 



z 1 V RT 

where 1/Z is the normalization factor. Let us select here the normalization Z — 1 and write: 



F = V/(c,r) + W2rJ^Yln( 



(50) 



We assume that the standard equilibrium concentrations ^*(c, T) are much smaller than the 
concentrations of A^. It is always possible because functions Uj are defined up to an additive constant. 

The formula for free energy is necessary to define the fast equilibria (45). Such an equilibrium is the 
minimizer of the free energy on the straight line parameterized by a: Ci = c° — auji, % = a. 

If we neglect the products qjdqUc, T)/dci as the second order small quantities then the minimizers 
have the very simple form: 

//i(c, T) 

#3 = 2^ U H p T ( 51 ) 



or 



where 



,;(c,T)exp( E -^ (c ' T) ) (52) 
df(c,T) 



is the chemical potential of Ai and 



= In 



{RTdj — is the chemical potential of -Bj). 

The thermodynamic equilibrium of the system of reactions Bj — > B t that corresponds to the 
reactions (46) is the free energy minimizer under given values of the conservation laws. 

For the systems with fixed volume, the stoichiometric conservation laws of the monomolecular system 
of reaction are sums of the concentrations of Bj which belong to the connected components of the 
reaction graph. Under the hypothesis of weak reversibility there is no other conservation law. Let the 
graph of reactions Bj — > B[ have d connected components C s and let V s be the set of indexes of those 
Bj which belong to C s : Bj E C s if and only if j £ V s . For each C s there exists a stoichiometric 
conservation law 

0s = = const (53) 

For any set of positive values of /3 S (s = 1, . . . , d) and given c, T there exists a unique conditional 
maximizer ^ cq of the free energy (50): For the compound Bj from the sth connected component (j £ V s ) 
this equilibrium concentration is 



c 1= R S*( C ' T ) 
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The positive values of concentrations ^ are the equilibrium concentrations (54) for some values of (3 S 
if and only if for any s = 1 , . . . , d and all j, I E V s 

&i = #i (55) 

(§j = ln^/d*)). This means that compounds are in equilibrium in every connected component C s the 
chemical potentials of compounds coincide in each component C s . The system of equations (55) together 
with the equilibrium conditions (52) constitute the equilibrium of the systems. All the equilibria form a 
linear subspace in the space with coordinates fa/ RT (i — 1, . . . , n) and (j = 1, . . . , q). 

In the expression for the free energy (50) we do not assume anything special about free energy of 
the mixture of Ai. The density of this free energy, /(c, T), may be an arbitrary smooth function (later, 
we will add the standard assumption about convexity of f(c, T) as a function of c). For the compounds 
Bi, we assume that they form a very small addition to the mixture of Ai, neglect all quadratic terms in 
concentrations of Bi and use the entropy of the perfect systems, p lap, for this small admixture. 

This approach results in the explicit expressions for the fast equilibria (52) and expression of the 
equilibrium compound concentrations through the values of the stoichiometric conservation laws (54). 

5.4. Markov Kinetics of Compounds 

For the kinetics of compounds transformations Bj — > B u the same hypothesis of the smallness of 
concentrations leads to the only reasonable assumption: The linear (monomolecular) kinetics with the 
rate constant K t j > 0. This "constant" is a function of c, T: f%(c, T). The order of indexes at k is inverse 
to the order of them in reaction: «y = K^j. 

The master equation for the concentration of Bj gives: 

It is necessary to find when this kinetics respect thermodynamics, i.e., when the free energy decreases due 
to the system (56). The necessary and sufficient condition for matching the kinetics and thermodynamics 
is: The standard equilibrium <;* (49) should be an equilibrium for (56), that is, for every j = 1, . . . , q 

= Kl ^*i (57) 

This condition is necessary because the standard equilibrium is the free energy minimizer for given c, T 
and J2j 9 = S*" ^ ne sum % conserves due to (56). Therefore, if we assume that F decreases 
monotonically due to (56) then the point of conditional minimum of F on the plane J2j S = const 
(under given c, T) should be an equilibrium point for this kinetic system. This condition is sufficient due 
to the Morimoto //-theorem (see Appendix 2). 

For a weakly reversible system, the set of the conditional minimizers of the free energy (54) coincides 
with with the set of positive equilibria for the master equations (56) (see Equation (132) in Appendix 2). 
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5. 5. Thermodynamics and Kinetics of the Extended System 

In this section, we consider the complete extended system, which consists of species Ai (i — 1, . . . , n) 
and compounds Bj (j = 1, . . . , q) and includes reaction of equilibration (45) and transformations of 
compounds Bj — > B\ which correspond to the reactions (46). 

Thermodynamic properties of the system are summarized in the free energy function (50). For kinetics 
of compounds we accept the Markov model (56) with the equilibrium condition (57), which guarantees 
matching between thermodynamics and kinetics. 

For the equilibration reactions (45) we select a very general form of the kinetic law. The only 
requirement is: This reaction should go to its equilibrium, which is described as the conditional 
minimizer of free energy F (52). For each reaction Qj ^ Bj (where the complex is a formal 
combination: Qj = £\ UjiAi) we introduce the reaction rate Wj. This rate should be positive if 

i 

and negative if 

RT 



^>5>T (»> 



The general way to satisfy these requirement is to select q continuous function of real variable Wj(x), 
which are negative if x > and positive if x < 0. For the equilibration rates we take 



Wj = Wj 



If several dynamical systems defined by equations x = Ji, ... x = J v on the same space have the 
same Lyapunov function F, then for any conic combination J = J2 k dkJk (a*; > 0, J2 k > 0) the 
dynamical system x = J also has the Lyapunov function F. 

The free energy (50) decreases monotonically due to any reaction Qj ^ Bj with reaction rate Wj (60) 
and also due to the Markov kinetics (56) with the equilibrium condition (57). Therefore, the free energy 
decreases monotonically due to the following kinetic system: 

dc 



at 

where the coefficients Kji satisfy the matching condition (57). 

This general system (61) describes kinetics of extended system and satisfies all the basic conditions 
(thermodynamics and smallness of compound concentrations). In the next sections we will study the QE 
approximations to this system and exclude the unknown functions Wj from it. 

5.6. QE Elimination of Compounds and the Complex Balance Condition 



In this section, we use the QE formalism developed for chemical kinetics in Section 4 for 
simplification of the compound kinetics. 
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First of all, let us describe L L , where the space L is the subspace in the extended concentration space 
spanned by the stoichiometric vectors of fast equilibration reactions (45). The stoichiometric vector for 
the equilibration reactions have a very special structure (47). Dimension of the space L is equal to the 
number of complexes: dimL = q. Therefore, dimension of L 1 - is equal to the number of components 
Af. dimL^ = n. For each Ai we will find a vector bi E L 1 - that has the following first n coordinates: 
bik — $ik for k — 1, . . . , n. The condition (pi, gj) = gives immediately: &i, n +j = Vj\. Finally, 



bi = ( 0, ■ ■ . , 0, 1 , 0, . . . , 0, V U , U 2i , • • • , Vqi) (62) 
i 

The corresponding slow variables are 

bi(c, q) = Cj + } j qjUji (63) 

In the QE approximation all Wj = and the kinetic equations (61) give in this approximation 

dbi v-^ 

In these equations, we have to use the dependence q(b). Here we use the QSS Michaelis and Menten 
assumption: The compounds are present in small amounts 

(k > q-j 

In this case, we can take bi instead of q (i.e., take fx(b,T) instead of /x(c, T)) in the formulas for 
equilibria (52): 

9 = S T ) ex P I ^ ) (65) 

In the final form of the QE kinetic equation there remain two "offprints" of the compound kinetics: 
Two sets of functions qj{b, T) > and Kji(b, T) > 0. These functions are connected by the identity (57). 
The final form of the equations is 

dt" = 1^ [ Kjl< " 1 ( ' > 6XP l RT J ~ ^ ' ^ 6XP V RT j j ji (66) 

The identity (57), Yli i^j K ji^i = J2i i^j P rov ides a sufficient condition for decreasing of free 
energy due to the kinetic equations (66). This is a direct consequence of two theorem: The theorem 
about the preservation of entropy production in the QE approximations (see Section 2 and Appendix 1) 
and the Morimoto //-theorem (see Appendix 2). Indeed, in the QE state the equilibrated reactions (45) 
Qj ^ Bj do not produce entropy and all changes in the total free energy are caused by the Markov 
kinetics Bi — > Bj. Due to the Morimoto //-theorem this change is negative: The Markov kinetics 
decrease the perfect free energy of compounds and do not affect the free energy of A4. In the QE 
approximation, the concentrations of Ai are changing together with concentrations of Bj because of 
the equilibrium conditions for reactions Qj ^ Bj. Due to the theorem of preservation of the entropy 
production, the time derivative of the total free energy in this QE dynamics coincides with the time 
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derivative of the free energy of Bj due to Markov kinetics. In addition to this proof, in Section 6 below 
we give the explicit formula for entropy production in (66) and direct proof of its positivity. 

Let us stress that the functions <;*(&, T) and Kji(b, T) participate in equations (66) and in identity (57) 
in the form of the product. Below we use for this product a special notation: 

<p jt (b, T) = Kjl (b, T)^(b, T) {j ^ I) (67) 

We call this function ipji(b, T) the kinetic factor. The identity (57) for the kinetic factor is 

<Pfl (6, T) = <Pij (b,T) for all j (68) 

We call the thermodynamic factor (or the Boltzmann factor) the second multiplier in the reaction rates 

iliip, T) = exp — (69) 



RT 

In this notation, the kinetic equations (66) have a simple form 

^ = (<Pji(b,T)^i(b,T) - ^(6,T)%(fo,r))^ (70) 

The general equations (70) have the form of "sum over complexes". Let us return to the more usual 
"sum over reactions" form. An elementary reaction corresponds to the pair of complexes 6;, Qj (46). 
It has the form 0/ — > Qj and the reaction rate is r = (pjifti. In the right hand side in (70) this reaction 
appears twice: first time with sign "+" and the vector coefficient Uj and the second time with sign "— " 
and the vector coefficient v\. The stoichiometric vector of this reaction is 7 = Vj — v%. Let us enumerate 
the elementary reactions by index p, which corresponds to the pair (j, I). Finally, we transform (46) into 
the sum over reactions form 



lj,lft 

= ^2vp(b,T)Q p {b,T)-f pi 
p 

In the vector form it looks as follows: 

db 



(71) 



dt 



Y,V P (b,T)n p (b,T) lp (72) 
p 

5. 7. The Big Michaelis-Menten-Stueckelberg Theorem 



Let us summarize the results of our analysis in one statement. 

Let us consider the reaction mechanism illustrated by Figure 2 (46): 

i i 

under the following asymptotic assumptions: 
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1. Concentrations of the compounds B p are close to their quasiequilibrium values (65) 

* = (1 + ^ = (1 + T) exp ( E< ^ (6 ' r) ) , * « 1 

(this may be due to the fast reversible reactions in (46)); 

2. Concentrations of the compounds B p are much smaller than the concentrations of the components 
Ai\ There is a small positive parameter £<1, q* = e£* and £* do not depend on s; 

3. Kinetics of transitions between compounds — >■ Bj is linear (Markov) kinetics with reaction rate 
constants kji = 

Under these assumptions, in the asymptotic 5, e — > 0, 5, e > kinetics of components At may be 
described by the reaction mechanism 

i i 

with the reaction rates 

where the kinetic factors cp p satisfy the condition (68): 

p,a p =v p, /3 p =v 

for any vector v from the set of all vectors {a p , f3 p }. This statement includes the generalized mass action 
law for r p and the balance identity for kinetic factors that is sufficient for the entropy growth as it is 
shown in the next Section 6. 



6. General Kinetics and Thermodynamics 

6. 1 . General Formalism 



To produce the general kinetic law and the complex balance conditions, we use "construction staging": 
The intermediate complexes with fast equilibria, the Markov kinetics and other important and interesting 
physical and chemical hypothesis. 

In this section, we delete these construction staging and start from the forms (69), (72) as the basic 
laws. We use also the complex balance conditions (68) as a hint for the general conditions which 
guarantee accordance between kinetics and thermodynamics. 

Let us consider a domain U in n-dimensional real vector space E with coordinates N±, . . . , N n . For 
each Ni a symbol (component) A { is given. A dimensionless entropy (or free entropy, for example, 
Massieu, Planck, or Massieu-Planck potential which correspond to the selected conditions [52]) S(N) 
is defined in U. "Dimensionless" means that we use S/R instead of physical S. This choice of units 
corresponds to the informational entropy (p \np instead of k-^p hip). 

The dual variables, potentials, are defined as the partial derivatives of S: 

* - ~m < 73 > 
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Warning: This definition differs from the chemical potentials (22) by the factor 1/RT: For constant 
volume the Massieu-Planck potential is — F/T and we, in addition, divide it on R. On the other hand, we 
keep the same sign as for the chemical potentials, and this differs from the standard Legendre transform 
for S. (It is the Legendre transform for function — S). 

The reaction mechanism is defined by the stoichiometric equations (17) 

i i 

(p = 1, . . . , m). In general, there is no need to assume that the stoichiometric coefficients a pi , (5 pi 
are integers. 

The assumption that they are nonnegative, a pi > 0, (3 P i > 0, may be needed to prove that the kinetic 
equations preserve positivity of iVj. If iVj is the number of particles then it is a natural assumption but 
we can use other extensive variables instead, for example, we included energy in the list of variables to 
describe the non-isothermal processes [53]. In this case, the coefficient au for the energy component 
A v in an exothermic reaction is negative. 

So, for variables that are positive (bounded from below) by their physical sense, we will use the 
inequalities a pi > 0, (3 P i > 0, when necessary, but in general, for arbitrary extensive variables, we do not 
assume positivity of stoichiometric coefficients. As it is usually, the stoichiometric vector of reaction is 
1 P = ftp — otp (the "gain minus loss" vector). 

For each reaction, a nonnegative quantity, reaction rate r p is defined. We assume that this quantity 
has the following structure: 

r p = y? p exp(a p ,/2) (74) 

where (a p ,p) = Ei^Ai- 

In the standard formalism of chemical kinetics the reaction rates are intensive variables and in kinetic 
equations for N an additional factor — the volume — appears. For heterogeneous systems, there may be 
several "volumes" (including interphase surfaces). 

Each reaction has it own "volume", an extensive variable V p (some of them usually coincide), and we 
can write 

= v p~tpvp ex p(«p> a) ( 75 ) 
p 

In these notations, both the kinetic and the Boltzmann factors are intensive (and local) characteristics 
of the system. 

Let us, for simplicity of notations, consider a system with one volume, V and write 

p 

Below we use the form (76). All our results will hold also for the multi-volume systems (75) under 
one important assumption: The elementary reaction 



i i 
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goes in the same volume as the reverse reaction 

/3 P iAi ->■ ^2 a Pi A i 

i i 

or symbolically 

V; = V- (77) 

If this condition (77) holds then the detailed balance conditions and the complex balance conditions will 
hold separately in all volumes V p . 

An important particular case of (76) gives us the Mass Action Law. Let us take the perfect free entropy 



s = -Y,NiM%)- 1 ) (78) 



Ci 

where q = Ni/V > are concentrations and c* > are the standard equilibrium concentrations. Under 
isochoric conditions, V = const, there is no difference between the choice of the main variables, N or c. 
For the perfect function (78) 

fa = ln , exp(a p ,/2) = Yl (79) 

and for the reaction rate function (74) we get 



The standard assumption for the Mass Action Law in physics and chemistry is that if and c* are functions 
of temperature: tp p = f p (T) and c* = c*(T). To return to the kinetic constants notation (20) we 
should write: 

Titer* p 

Equation (76) is the general form of the kinetic equation we would like to study. In many senses, this 
form is too general before we impose restrictions on the values of the kinetic factors. For physical and 
chemical systems, thermodynamics is a source of restrictions: 

1 . The energy of the Universe is constant. 

2. The entropy of the Universe tends to a maximum. 
(R. Clausius, 1865 [54].) 

The first sentence should be extended: The kinetic equations should respect several conservation laws: 
Energy, amount of atoms of each kind (if there is no nuclear reactions in the system) conservation of total 
probability and, sometimes, some other conservation laws. All of them have the form of conservation 
of values of some linear functionals: l(N) = const. If the input and output flows are added to the 
system then 

H/fAfl . . 

: Vv m l m - v out l(N) 



dt 
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where v m ' out are the input and output fluxes per unit volume, l m are the input densities (concentration). 
The standard requirement is that every reaction respects all these conservation laws. The formal 
expression of this requirement is: 

l{j p ) = for all p (81) 

There is a special term for this conservation laws: The stoichiometric conservation laws. All the main 
conservation laws are assumed to be the stoichiometric ones. 

Analysis of the stoichiometric conservation laws is a simple linear algebra task: We have to find 
the linear functionals that annulate all the stoichiometric vectors j p . In contrast, entropy is not a linear 
function of N and analysis of entropy production is not so simple. 

In the next subsection we discuss various conditions which guarantee the positivity of entropy 
production in kinetic equations (76). 

6.2. Accordance Between Kinetics and Thermodynamics 

6.2.1. General Entropy Production Formula 

Let us calculate dS/dt due to equations (76): 

dS _ v ds dNi 

~dt ~ ^ dNi dt 

i 

= - ^ piV IpiVp exp(a p , p) (82 ) 

i p 

= -V5^(7p,/x)^ p exp(a p ,/2) 
p 

An auxiliary function 9(X) of one variable A E [0, 1] is convenient for analysis of dS/dt (it was 
studied by Rozonoer and Orlov [55], see also [25]: 

0(A) = ]T <p p exp[(fi, (Xa p + (1 - \)B p ))) (83) 
p 

With this function, the entropy production (82) has a very simple form: 



dS__ y d0(A) 



dt dA 



(84) 

A=l 



The auxiliary function 9(\) allows the following interpretation. Let us introduce the deformed 
stoichiometric mechanism with the stoichiometric vectors, 



a p (X) = \a p + (1 - \)B p , p (X) = \/3 p + (1 - A) 



a, 



, which is the initial mechanism when A = 1, the inverted mechanism with interchange of a and (3 when 
A = 0, the trivial mechanism (the left and right hand sides of the stoichiometric equations coincide) 
when A = 1/2. 

For the deformed mechanism, let us take the same kinetic factors and calculate the Boltzmann factors 
with a p (\): 

r pW = ^ P exp(a p (A),/i) 
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In this notation, the auxiliary function 0(A) is a sum of reaction rates for the deformed 
reaction mechanism: 

0(A) = J> p (A) 
p 

In particular, 0(1) = J2 P r p> m i s i s j ust me sum of reaction rates. 
Function 0(A) is convex. Indeed 

= E Vp(7p, A) 2 exp[(/2, (\a p + (1 - A)/? p ))] > 
p 

This convexity gives the following necessary and sufficient condition for positivity of 
entropy production: 

— > if and only if 0(A) < 0(1) for some A < 1 

In several next subsections we study various important particular sufficient conditions for positivity 
of entropy production. 

6.2.2. Detailed Balance 

The most celebrated condition which gives the positivity of entropy production is the principle of 
detailed balance. Boltzmann used this principle to prove his famous //-theorem [27]. 
Let us join elementary reactions in pairs: 

MpiAi ^ Yl (85) 

i i 

After this joining, the total amount of stoichiometric equations decreases. If there is no reverse reaction 
then we can add it formally, with zero kinetic factor. We will distinguish the reaction rates and kinetic 
factors for direct and inverse reactions by the upper plus or minus: 

r t = ft exp(a p , p) , fZ = ip~ exp(/3 p , p) , r p = r+ - r~ 



p 

In this notation, the principle of detailed balance is very simple: The thermodynamic equilibrium 
in the direction j p , given by the standard condition (j p , jx) = 0, is equilibrium for the corresponding 
pair of mutually reverse reactions from (85). For kinetic factors this transforms into the simple and 
beautiful condition: 

(ft exp(a p , fi) = (f~ exp(/3 p , fi) & (7 P , fi) = 

therefore 

ft = f^ (8V) 
For the systems with detailed balance we can take y? p = (p p = (p~ and write for the reaction rate: 



r p = v9 p (exp(a p , p) - exp(/3 p , p)) 
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M. Feinberg called this kinetic law the "Marselin-De Donder" kinetics [46]. This representation of the 
reaction rates gives for the auxiliary function 0(A): 

0(A) = <Pp(exp[(fi, (Aap + (1 - \)P P ))] + exp[(/2, (A/3, + (1 - X)a p ))]) (88) 
p 

Each term in this sum is symmetric with respect to change A i — >■ (1 — A). Therefore, 0(1) = 0(0) and, 
because of convexity of 0(A), 0'(\) > 0. This means positivity of entropy production. 

The principle of detailed balance is a sufficient but not a necessary condition of the positivity of 
entropy production. This was clearly explained, for example, by L. Onsager [56,57]. Interrelations 
between positivity of entropy production, Onsager reciprocal relations and detailed balance were 
analyzed in detail by N.G. van Kampen [58]. 

6.2.3. Complex Balance 

The principle of detailed balance gives us 0(1) = 0(0) and this equality holds for each pair of mutually 
reverse reactions. 

Let us start now from the equality 0(1) = 0(0). We return to the initial stoichiometric equations (17) 
without joining the direct and reverse reactions. The equality reads 

^2 V P exp {ft , a p ) = ^2 Vp ex P (fi'Pp) ( 89 ) 
p p 

Exponential functions exp(/i, y) form linearly independent family in the space of functions of ft for 
any finite set of pairwise different vectors y. Therefore, the following approach is natural: Let us 
equalize in (89) the terms with the same Boltzmann-type factor exp(/i, y). Here we have to return to 
the complex-based representation of reactions (see Section 5.1). 

Let us consider the family of vectors {a p ,(3 p } (p = 1, . . . , m). Usually, some of these vectors 
coincide. Assume that there are q different vectors among them. Let yx, . . . , y q be these vectors. For 
each j — 1 , . . . , q we take 

Rj = {p I a p = Vj} , R j = {p\P P = Vj} 
We can rewrite the equality (89) in the form 



3=1 



P dR+ p€R~ 



(90) 



The Boltzmann factors exp(/t, yj) form the linearly independent set. Therefore the natural way to meet 
these condition is: For any j = l,...,q 

E v P - E vp = (91) 

per} peRj 

This is the general complex balance condition. This condition is sufficient for entropy growth, because 
it provides the equality 0(1) = 0(0). 
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If we assume that ip p are constants or, for chemical kinetics, depend only on temperature, then the 
conditions (91) give the general solution to equation (90). 

The complex balance condition is more general than the detailed balance. Indeed, this is obvious: 
For the master equation (56) the complex balance condition is trivially valid for all admissible constants. 
The first order kinetics always satisfies the complex balance conditions. On the contrary, the class of 
the master equations with detailed balance is rather special. The dimension of the class of all master 
equations has dimension n 2 — n (constants for all transitions Aj — > Aj are independent). For the 
time-reversible Markov chains (the master equations with detailed balance) there is only n(n + l)/2 — 1 
independent constants: n — 1 for equilibrium state and n(n — l)/2 for transitions — > A, (i > j), 
because for reverse transitions the constant can be calculated through the detailed balance. 

In general, for nonlinear reaction systems, the complex balance condition is not necessary for entropy 
growth. In the next section we will give a more general condition and demonstrate that there are systems 
that violate the complex balance condition but satisfy this more general inequality. 

6.2.4. G-Inequality 

Gorban [25] proposed the following inequality for analysis of accordance between thermodynamics 
and kinetics: 0(1) > 9(0). This means that for any values of fi 



^(^ p exp(/2,a p ) > ^^ p exp(/i,/3 p 



(92) 



p p 
In the form of sum over complexes (similarly to (90)) it has the form 



Y^p-Y^p 



peRI 



> 



(93) 



Let us call these inequalities, (92), (93), the G-inequalities. 

Here, two remarks are needed. First, functions exp(/2, yj) are linearly independent but this does not 
allow us to transform inequalities (93) similarly to (91) even for constant kinetic factors: Inequality 
between linear combinations of independent functions may exist and the "simplified system" 

Y V P ~ Y V P > for all j 
P eRj pen- 
is not equivalent to the G-inequality. 

Second, this simplified inequality is equivalent to the complex balance condition (with equality 
instead of >). Indeed, for any p — 1, . . . , m there exist exactly one ji and one j 2 ^ j\ with properties: 
p G Rj^, p E Rj 2 . Therefore, for any reaction mechanism with reaction rates (74) the identity holds: 



E 



v P - Yl 



pern 



If all terms in this sum are non-negative then all of them are zeros. 
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Nevertheless, if at least one of the vectors yj is a convex combination of others, 

^2 = yj for some Afc - °> Afc = 1 

fc, k^j k, k^j 

then the G inequality has more solutions than the condition of complex balance. Let us take a very 
simple example with two components, Ax and A 2 , three reactions and three complexes: 

2A X ^A Y + A 2 , 2A 2 ^ A x + A 2 , 1A X ^ 2A 2 

yi = (2,0),j/ 2 = (0,2),2/ 3 = (l,l), 
i?+ = {1,3}, i?+ = {2, -3}, i? 3 + = {-1, -2} 

i?r = {"I, -3}, R 2 = {-2, 3}, ^ = {1,2} 
The complex balance condition for this system is 

(V?l - (f-l) + (<f 3 - <Ps) = 
(<f 2 - y?_ 2 ) - (v? 3 - ^-3) = 

The G-inequality for this system is 

(</?! + tp z - v?_i - if^)a 2 + {ip 2 + v?_ 3 - y?_ 2 - V5 3 )& 2 
+ + tp-2 — fi — ^2)0^ > for all a,b > 

(here, a, b stand for exp(fii), exp(/2 2 )). Let us use for the coefficients at a 2 and b 2 notations ip a and 
ip b - Coefficient at ab in (95) is — (ip a + ip b ), linear combinations ^> a = tpx + ip 3 — tp_\ — <^_ 3 and 
4>b — V?2 + <^-3 — V-2 — ¥3 are linearly independent functions of variables (i = ±1, ±2, ±3) and we 
get the following task: To find all pairs of numbers (ip a , ipb) £ which satisfy the inequality 

ip a a 2 + ipbb 2 > (ip a + ipb)ab for all a, 6 > 

Asymptotics a — > and & — )■ give ^ > 0. 

Let us use homogeneity of functions in (95), exclude one normalization factor from a, b and one 
factor from ip a , tp b and reduce the number of variables: b — 1 — a, ip a — 1 — ipb- We have to find all such 
ipb e [0, 1] that for all a e]0, 1[ 

a 2 (l - + (1 - a)V& - o(l - 0) > 

The minimizer of this quadratic function of a is a min = \ + lip b , a min e]0, 1[ for all ipb G [0, 1]. 
The minimal value is — 2(|'0& — |) 2 . It is nonnegative if and only if ipb = |. When we return to the 
non-normalized variables ip a , ipb then we get the general solution of the G-inequality for this example: 
ipa = ipb > 0. For the kinetic factors this means: 

(<pi - + 2(v? 3 - v?_ 3 ) - (y? 2 - ^-2) = 

(pi - p_i) + (<P3 - tp-s) > (96) 

(<p 2 - ¥-2) - (<£>3 - <P-3) > 



(94) 



(95) 
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These conditions are wider (weaker) than the complex balance conditions for this example (94). 

In the Stueckelberg language [26], the microscopic reasons for the G-inequality instead of the 
complex balance (68) can be explained as follows: Some channels of the scattering are unknown 
(hidden), hence, instead of unitarity of ^-matrix (conservation of the microscopic probability) we have 
an inequality (the microscopic probability does not increase). 

We can use other values of A G [0, 1[ in inequality 0(1) > 9(\q) and produce constructive sufficient 
conditions of accordance between thermodynamics and kinetics. For example, condition 9(1) > 9(1/2) 
is weaker than 9(1) > 9(0) because of convexity 9(X). 

One can ask a reasonable question: Why we do not use directly positivity of entropy production 
(9'(1) > 0) instead of this variety of sufficient conditions. Of course, this is possible, but inequalities 
like 9(1) > 9(0) or equations like 9(1) = 9(0) include linear combinations of exponents of linear 
functions and often can be transformed in algebraic equations or inequalities like in the example above. 
Inequality ^'(l) > includes transcendent functions like /exp / (where / is a linear function) which 
makes its study more difficult. 

7. Linear Deformation of Entropy 

7.1. If Kinetics Does not Respect Thermodynamics then Deformation of Entropy May Help 

Kinetic equations in the general form (76) are very general, indeed. They can be used for the 
approximation of any continuous time dynamical system on compact U [59]. In previous sections we 
demonstrated how to construct the system in the form (76) with positivity of the entropy production 
when the entropy function is given. 

Let us consider a reverse problem. Assume that a system in the form (76) is given but the entropy 
production is not always positive. How to find a new entropy function for this system to guarantee the 
positivity of entropy production? 

Existence of such an entropy is very useful for analysis of stability of the system. For example, let us 
take an arbitrary Mass Action Law system (80). This is a rather general system with the polynomial right 
hand side. Its stability or instability is not obvious a priori. It is necessary to check whether bifurcations 
of steady states, oscillations and other interesting effects of dynamics are possible for this system. 

With the positivity of entropy productions these questions are much simpler (for application of 
thermodynamic potentials to stability analysis see, for example, [33,59-61]). If dS/dt > and it 
is zero only in steady states, then any motion in compact U converges to a steady state and all the 
non- wandering points are steady states. (A non- wandering point is such a point x E U that for any 
T > and e > there exists such a motion c(t) G U that (i) ||c(0) — x\\ < e and (ii) \\c(T') — x\\ < e 
for some T' > T: A motion returns in an arbitrarily small vicinity of x after an arbitrarily long time.) 
Moreover, the global maximizer of S in U is an asymptotically stable steady state (at least, locally). It is 
a globally asymptotically stable point if there is no other steady state in U . 

For the global analysis of an arbitrary system of differential equations, it is desirable either to construct 
a general Lyapunov function or to prove that it does not exist. For the Lyapunov functions of the general 
form this task may be quite difficult. Therefore, various finite-dimensional spaces of trial functions are 
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often in use. For example, quadratic polynomials of several variables provide a very popular class of 
trial Lyapunov function. 

In this section, we discuss the n-parametric families of Lyapunov functions which are produced by 
the addition of linear function to the entropy: 

S(N) ^ S Ap (N) = S(N) - A & N i ( 97 ) 

i 

The change in potentials p is simply the addition of Ap: pi \-± pi + A/V 

Let us take a general kinetic equation (76). We are looking for a transformation that does not change 
the reaction rates. The Boltzmann factor Vt p = exp(/2, a p ) transforms due to the change of the entropy: 
Vt p i — v Vt p exp(A/2, a p ). Therefore, to preserve the reaction rate, the transformation of the kinetic factors 
should be ip p H- ip p exp (A/2, a p ) in order to keep the product r p = Q p (f p constant. 

For the new entropy, S = Sa/i, with the new potential and kinetic factors, the entropy production is 
given by (98): 

p 

= ~ E(7p> ^ + A ^X d ex P(«P' ^ ld ) (98) 
p 

dS old 



^( 7p ,A/i)^ ld ex P (a p , / i old ) 



dt 

p 

where the superscript "old" corresponds to the non-deformed quantities. 
7. 2. Entropy Deformation for Restoration of Detailed Balance 



It may be very useful to find such a vector A/2 that in new variables <p+ = ip . For the analysis of 
the detailed balance condition, we group reactions in pairs of mutually inverse reactions (85). Let us 

"p = r p ~ V ¥p 



consider an equation of the general form (86) with r p = r+ — r , > 0. 



The problem is: To find such a vector Ap, that 

V?+ exp (Ap, a p ) = ip~ exp (Ap, f} p ) (99) 
or, in the equivalent form of the linear equation 

(Ap, lp )=\n{^j (100) 

The necessary and sufficient conditions for the existence of such Ap are known from linear algebra: 
For every set of numbers a p (p = 1, . . . , m) 

Y^a plp = 0^J2 a ? ln (j^) =0 ( 101 ) 
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To check these conditions, it is sufficient to find a basis of solutions of the uniform systems of 
linear equations 

J^a P 7pi = (i = 1, . . .,m) 
p 

(that is, to find a basis of the left kernel of the matrix T, coimr, where T = (7^)) and then check for 
these basis vectors the condition J2 P a p m (i^r~) = to prove or disprove that the vector with coordinates 

In belongs to the image of T, imT. 

For some of the reaction mechanisms it is possible to restore the detailed balance condition for the 
general kinetic equation unconditionally. For these reactions, for any set of positive kinetic factors, 
there exists such a vector A/2 that the detailed balance condition (100) is valid for the deformed entropy. 
According to (101) this means that there is no nonzero solution a p for the equation ^ p a p 7 p = 0. In 
other words, vectors 7 P are independent. 

7. 3. Entropy Deformation for Restoration of Complex Balance 

The complex balance conditions (91) are, in general, weaker than the detailed balance but they are 
still sufficient for the entropy growth. 

Let us consider an equation of the general form (86). We need to find such a vector A/2 
that in new variables with the new entropy and kinetic factors the complex balance conditions 

Z P eRj - E P6RJ ^ eW = hold " 

For our purpose, it is convenient to return to the presentation of reactions as transitions between 

complexes. The complexes, 9 1; . . . , Q q are the linear combinations, Qj = (yj, A). 

Each elementary reaction (17) with the reaction number p may be represented in the form Oj — > 6;, 

where Qj = UjAj, p G Rj (a p = yf) and p G Rj (/3 P = y{). For this reaction, let us use the notation 

ip p = <pij. We used this notation in the analysis of kinetics of compounds (Section 5.6). The complex 

balance conditions are 

^2{<Pi3-<Pji) = (102) 
To obtain these conditions after the entropy deformation, we have to find such A/i that 

^2 fan ex P ( A A Vi) ~ fji ex P Vi)) = ( 103 ) 

This is exactly the equation for equilibrium of a Markov chain with transition coefficients (fij. 
Vector (A/i, yj) should be an equilibrium state for this chain (without normalization to the unit sum 
of coordinates). 

For this finite Markov chain a graph representation is useful: Vertices are complexes and oriented 
edges are reactions. To provide the existence of a positive equilibrium we assume weak reversibility of 
the chain: If there exists an oriented path from Oj to 6; then there exists an oriented path from 6; to Qj. 

Let us demonstrate how to transform this problem of entropy deformation into a linear algebra 
problem. First of all, let us find any positive equilibrium of the chain, q* > 0: 

- w?) = o (104) 
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This is a system of linear equations. If we have already an arbitrary equilibrium of the chain then 
other equilibria allow a very simple description. We already found this description for kinetics of 
compounds (54) 

Let us consider the master equation for the Markov chain with coefficients cp t j and apply the formalism 
from Appendix 2: 

Let the graph of complex transformations Qj — > Qi have d connected components C s and let V s be 
the set of indexes of those Qj which belong to C s : Qj 6 C s if and only if j E V s . For each C s there 
exists a conservation law for the master equation (53), /? s (<r) = J2jev s 9- 

For any set of positive values of j3 a (s = l,...,q) there exists a unique equilibrium vector <^ eq 
for (105) with this values (3 S (54), (132). The set of equilibria is a linear space with the natural coordinates 
/3 S (s = 1, . . . , d). We are interested in the positive orthant of this space, f3 s > 0. For positive /3 S , 
logarithms of <^ eq form a (/-dimensional linear manifold in R q (133). The natural coordinates on this 
manifold are In (3 S . 

Let us notice that the vector ^° with coordinates 

is also an equilibrium for (105). This equilibrium is normalized to unit values of all S (^°). In the 
coordinates In (3 S this is the origin. The equations for A/2 are 

(A/2, V j) - hxfis = In S ° for j e V s (106) 

This is a system of linear equations with respect to n + d variables A/ij (i = l,...,n) and 
In (3 S (s = 1, . . . , d). Let the coefficient matrix of this system be denoted by M. 

Analysis of solutions and solvability of such equations is one of the standard linear algebra tasks. If 
this system has a solution then the complex balance in the original system can be restored by the linear 
deformation of the entropy. If this system is solvable for any right hand side, then for this reaction 
mechanism we always can find the entropy, which provides the complex balance condition. 

Unconditional solvability of (106) means that the left hand side matrix of this system has rank q. Let 
us express this rank through two important characteristics: It is rank{7i, . . . , 7 m } + d, where d is the 
number of connected components in the graph of transformation of complexes. 

To prove this formula, let us write down the matrix M of the system (106). First, we change the 
enumeration of complexes. We group the complexes from the same connected component together and 
arrange these groups in the order of the connected component number. After this change of enumeration, 
{1, . . . , 1^1} = V u {\Vx\ + 1, . . ■ , |Vi| + \V 2 \} = V 2 , {\V 1 \ + \V 2 \ + ... + +1,...,\V 1 \ + \V 2 \ + 
■ ■■ + \V d \} = V d . 
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Let yj be here the row vector. The matrix is 



M 



Vi 

y\vi\ 

. 2/|Vi|+... 



1 ... 



-I Veil 



1 
1 










(107) 



M consists of d blocks M s , which correspond to connected components C s of the graph of 
transformation of complexes: 



2/|V 1 | + ...+|y s _ 1 |+i 



y\v x \+...+\Vs\ o 



(108) 



The first n columns in this matrix are filled by the vectors yj of complexes, which belong to the 
component C s , then follow s — 1 columns of zeros, after that, there is one column of units, and then 
again zeros. Here, in (107), (108) we multiplied the last d columns by —1. This operation does not 
change the rank of the matrix. 

Other elementary operations that do not change the rank are: We can add to any row (column) a linear 
combination of other rows (columns). 

We will use these operations to simplify blocks (108) but first we have to recall several properties of 
spanning trees [62] . Let us consider a connected, undirected graph G with the set of vertices V and the 
set of edges £ C V x V. A spanning tree of G is a selection of edges of G that form a tree spanning every 
vertex. For a connected graph with V vertices, any spanning tree has V — 1 edges. Let for each vertex 
Qj of G a n-dimensional vector y t is given. Then for every edge (Qj, Qi) G £ a vector 7^ = yj — y x 
is defined. We identify vectors 7 and —7 and the order of j, I is not important. Let us use r G for this 
set of 7j 7 : 

T G = {y J -y l \(Q J ,Q l )e£} 
For any spanning tree T of graph G we have the following property: 



spanr G = spanlY 



(109) 



in particular, rankr G = ranklY. 

For the digraphs of reactions between complexes, we create undirected graphs just by neglecting the 
directions of edges. We keep for them the same notations as for original digraphs. Let us select any 
spanning tree T s for the connected component C s in the graph of transformation of complexes. In T s we 
select arbitrarily a root complex. After that, any other complex Qj in C s has a unique parent. This is the 
vertex connected to it on the path to the root. For the root complex of C s we use special notation 9°. 
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Now, we transform the block (108) without change of rank: For each non-root complex we subtract 
from the corresponding row the row which correspond to its unique parent. After these transformations 
(and, may be, some permutations of rows), the block M s get the following form: 

' it o 

7i 

. y° s o 

Here, {7^, 7|, . . . , 7j s v - 1 |_ 1 } is T Ts for the spanning tree T s and y° is the coefficient vector for the root 
complex 9°. 

From the obtained structure of blocks we immediately find that the rank of the rows with 7 is 
rank{7i, . . . ,7 m } + d due to (109). Additional d rows with y° are independent due to their last 
coordinates and add d to rank. Finally 

rankM = rank{7i, . . . ,7 m } + d (111) 

Obviously, rankM < q. 

In particular, from the formula (1 1 1) immediately follows the description of the reaction mechanisms, 
for which it is always possible to restore the thermodynamic properties by the linear deformation of the 
entropy. 

The deficiency zero theorem. If rankM = q then it is always possible to restore the positivity of the 
entropy production by the linear deformation of the entropy. 

Feinberg [63] called the difference q — rankM the deficiency of the reaction network. For example, 
for the "Michaelis-Menten" reaction mechanism E + S ^ ES =^ P + S rank{7i,7 2 } = 2, d — 1, 
q — 3, rankM = 3 and deficiency is 0. 

For the adsorption (the Langmuir-Hinshelwood) mechanism of CO oxidation (38) 
rank{7i, 72, 73} = 3, d = 3, q = 6, rankM = 6 and deficiency is 0. To apply the results 
about the entropy deformation to this reaction mechanism, it is necessary to introduce an inverse 
reaction to the third elementary reaction in (38), PtO+PtCO— 7-C02+2Pt with an arbitrarily small but 
positive constant in order to make the mechanism weakly reversible. 

Let us consider the Langmuir-Hinshelwood mechanism for reduced list of components. Let us assume 
that the gas concentrations are constant because of control or time separation or just as a model "fast" 
system and just include them in the reaction rate constants for intermediates. Then the mechanism 
is 2Pt=^2PtO, Pt=±PtCO, PtO+PtC0^2Pt. For this system, rank{7i, 72, 73} = 2, d = 2, q = 5, 
rankM = 4 and deficiency is 1. Bifurcations in this system are known [33]. 

For the fragment of the reaction mechanism of the hydrogen combustion (44), rank{7i, . . . , 7 m } = 6, 
d = 7, q = 16, rankM = 6 + 7=13 and deficiency is 3. 








1 



(110) 
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7.4. Existence of Points of Detailed and Complex Balance 

Our formulation of the conditions of detailed and complex balance is not standard: We formulate 
them as the identities (87) and (89). These identities have a global nature and describe the properties of 
reaction rates for all states. 

The usual approach to the principle of detailed balance is based on equilibria. The standard 
formulation is: In all equilibria every process is balanced with its reverse process. Without special 
forms of kinetic law this principle cannot have any consequences for global dynamics. This is a trivial 
but not widely known fact. Indeed, let a system c = F(c) be given in a domain U cM. n and 71, . . . , j n is 
an arbitrary basis in W l . In this basis, we can always write: F(c) = YT P =i r p( c )lp- F° r an Y equilibrium 
c*, r p (c*) = 0. All the "reaction rates" r p (c) vanish simultaneously. This "detailed balance" means 
nothing for dynamics because F is an arbitrary vector field. Of course, if the system of vectors {7^} is 
not a basis but any complete system of vectors then such "detailed balance" conditions, r p (c*) = 0, also 
do not imply any specific features of dynamics without special hypotheses about functions r p (c). 

Nevertheless, if we fix the kinetic law then the consequences may be very important. For example, if 
kinetics of elementary reactions follow the Mass Action law then the existence of a positive equilibrium 
with detailed balance implies existence of the Lyapunov function in the form of the perfect free entropy: 



where c* is that positive equilibrium with detailed balance (see, for example, [33]). 

In this section we demonstrate that for the general kinetic law (74), which gives the expression 
of reaction rates through the entropy gradient, if the kinetic factors are constant (or a function of 
temperature) then the existence of the points of detailed (or complex) balance means that the linear 
deformation of the entropy exists which restores the global detailed (or complex) balance conditions (87) 



The condition that the kinetic factors are constant means that for a given set of values {ip p } a state 
with any admissible values of fx is physically possible (admissible). This condition allows us to vary the 
potentials fi independently of {(f p }. 

Let us assume that for the general kinetic system with the elementary reaction rates given by (74) a 
point of detailed balance exists. This means that for some value of fx — fx* (the detailed balance point in 
the Legendre transform) and for all p r+ = r~: 



This formula is exactly the condition (99) of existence of A/2 which allow us to deform the entropy for 
restoring the detailed balance in the global form (87). 

If we assume that the point of complex balance exists then there exists such a value of fx — fx* (a point 
of complex balance in the Legendre transform) that 




(or (89)). 



V9+exp(a p ,/2* 



ip p exp(a p ,/2* 




This is exactly the deformation condition (103) with A/2 = fx*. 
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To prove these statements we used an additional condition about possibility to vary fi under 
given {cp p }. 

So, we demonstrated that for the general kinetic law (74) the existence of a point of detailed balance 
is equivalent to the existence of such linear deformation of the entropy that the global condition (87) 
holds. Analogously, the existence of a point of complex balance is equivalent to the global condition of 
complex balance after some linear deformation of the entropy. 

7.5. The Detailed Balance is Needed More Often than the Complex Balance 

The complex balance conditions are mathematically nice and more general than the principle of 
the detailed balance. They are linked by Stueckelberg to the Markov models ("S'-matrix models") of 
microscopic kinetics. Many systems satisfy these conditions (after linear deformation of the entropy) 
just because of the algebraic structure of the reaction mechanism (see Section 7.3). Nevertheless, it is 
used much less than the classical detailed balance. 

The reason for the rare use of complex balance is simple: It is less popular because the stronger 
condition, the principle of detailed balance, is valid for most of physical and chemical systems. Onsager 
revealed the physical reason for detailed balance [56,57]. This is microreversibility: The microscopic 
laws of motion are invertible in time: If we observe the microscopic dynamics of particles in the 
backward movie then we cannot find the difference from the real world. This difference occurs in the 
macroscopic world. 

In microphysics and the S'-matrix theory this microreversibility property has the name "T-invariance". 

Let us demonstrate how T-invariance in micro-world implies detailed balance in macro-world. 

Following Gibbs, we accept the ensemble-based point of view on the macroscopic states: They are 
probability distributions in the space of detailed microscopic states. 

First of all, we assume that under given values of conservation laws equilibrium state exists and 
is unique. 

Second assumption is that the rates of elementary processes are microscopically observable quantities. 
This means that somebody (a "demon"), who observes all the events in the microscopical world can count 
the rates of elementary reactions. 

Because of T-invariance and uniqueness of equilibrium, the equilibrium is T-invariant: If we change 
all the microscopic time derivatives (velocities) v to — v then nothing will change. 

T-transformation changes all reactions to the reverse reactions, just by reversion of arrows, but the 
number of the events remains the same: Any reaction transforms into its reverse reaction but does not 
change the reaction rate. This can be formulated also as follows: T-transformation maps all r+ into the 
corresponding r~. 

Hence, because of the T-invariance, the equilibrium rate of each reaction is equal to the equilibrium 
rate of the reverse reaction. 

The violation of uniqueness of equilibrium for given values of conservation laws seems improbable. 
Existence of several equilibria in thermodynamics is quite unexpected for homogeneous systems 
but requires more attention for the systems with phase separation. Nevertheless, if we assume 
that a multi-phase system consists of several homogeneous phases, and each of these phases is 
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in uniform equilibrium, then we return to the previous assumption (with some white spots for 
non- uniform interfaces). 

T-invariance may be violated if the microscopic description is not reversible in time. Magnetic 
field and the Coriolis force are the classical examples for violation of the microscopic reversibility. 
In a linear approximation near equilibrium the corresponding modification of the Onsager relations 
give the Onsager-Casimir relations [64]. There are several attempts for nonlinear formulation of the 
Onsager-Casimir relations (see [65]). 

The principle of detailed balance seems to be still the best nonlinear version of the Onsager relations 
for T-invariant systems, and the conditions of the complex balance seem to give the proper relations 
between kinetic coefficients in the absence of the microscopic reversibility for nonlinear systems. It is 
important to mention here that all these relations are used together with the general kinetic law (74). 

Observability of the rates of elementary reactions deserves a special study. Two approaches to the 
reaction rate are possible. If we accept that the general kinetic law (74) is valid then we can find the 
kinetic factors by observation of dc/dt in several points because the Boltzmann factors are linearly 
independent. In this sense, they are observable but one can claim the approximation point of view and 
state that the general kinetic law (74) without additional conditions on kinetic factors is very general 
and allows to approximate any dynamical system. From this point of view, kinetic coefficients are 
just some numbers in the approximation algorithm and are not observable. This means that there is 
no such a microscopic thing as the rate of elementary reaction, and the set of reactions serves just for 
the approximation of the right hand side of the kinetic equation. We cannot fully disprove this point 
of view but can just say that in some cases the collision-based approach with physically distinguished 
elementary reactions is based on the solid experimental and theoretical background. If the elementary 
reactions physically exist then the detailed balance for T-invariant systems is proved. 

8. Conclusions 

We present the general formalism of the Quasiequilibrium approximation (QE) with the proof of the 
persistence of entropy production in the QE approximation (Section 2). 

We demonstrate how to apply this formalism to chemical kinetics and give several examples for 
the Mass Action law kinetic equation. We discuss the difference between QE and Quasi-Steady-State 
(QSS) approximations and analyze the classical Michaelis-Menten and Briggs-Haldane model reduction 
approaches (Section 3). After that, we use ideas of Michaelis, Menten and Stueckelberg to create a 
general approach to kinetics. 

Let us summarize the main results of our discussion. First of all, we believe that this is the finish 
of the Michaelis-Menten-Stueckelberg program. The approach to modeling of the reaction kinetics 
proposed by Michaelis and Menten in 1913 [8] for enzyme reactions was independently in 1952 applied 
by Stueckelberg [26] to the Boltzmann equation. 

The idea of the complex balance (cyclic balance) relations was proposed by Boltzmann as an answer 
to the Lorentz objections against Boltzmann's proof of the //-theorem. Lorentz stated that the collisions 
of the polyatomic molecules may have no inverse collisions. Cercignani and Lampis [28] demonstrated 
that the Boltzmann //-theorem based on the detailed balance conditions is valid for the polyatomic 
molecules under the microreversibility conditions and this new Boltzmann's idea was not needed. 
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Nevertheless, this seminal idea was studied further by many authors [29-31] mostly for linear systems. 
Stueckelberg [26] proved these conditions for the Boltzmann equation. He used in his proof the S-matrix 
representation of the micro-kinetics. 

Some consequences of the Stueckelberg approach were rediscovered for the Mass Action law kinetics 
by Horn and Jackson in 1972 [51] and supplemented by the "zero deficiency theorem" [63]. This is 
the history. 

In our work, we develop the Michaelis-Menten-Stueckelberg approach to general kinetics. This 
is a combination of the QE (fast equilibria) and the QSS (small amounts) approaches to the real or 
hypothetical intermediate states. These intermediate states (compounds) are included in all elementary 
reactions (46) as it is illustrated in Figure 1. Because of the small amount, the free energy for these 
compounds Bi is perfect (48), the kinetics of compounds is the first order Markov kinetics and satisfies 
the master equation. 

After that, we use the combination of QE and QSS approximations and exclude the concentrations 
of compounds. For the general kinetics the main result of this approach is the general kinetic law (74). 
Earlier, we just postulated this law because of its convenient and natural form [25,53], now we have the 
physical framework where this law can be proved. 

We do not assume anything about reaction rates of the main reactions (17). We use only 
thermodynamic equilibrium, the hypothesis about fast equilibrium with compounds and the smallness 
of concentration of compounds. This smallness implies the perfect entropy and the first order kinetics 
for compounds. After that, we get the reaction rate functions from the qualitative assumptions about 
compounds and the equilibrium thermodynamic data. 

For example, if we relax the assumption about fast equilibrium and use just smallness of compound 
concentrations (the Briggs-Haldane QSS approach [12-14]) then we immediately need the formulas for 
reaction rates of compound production. Equilibrium data become insufficient. If we relax the assumption 
about smallness of concentrations then we lose the perfect entropy and the first order Markov kinetics. 
So, only the combination of QE and QSS gives the desired result. 

For the kinetics of rarefied gases the mass action law for elastic collisions (the Boltzmann equation) 
or for inelastic processes like chemical reactions follows from the "molecular chaos" hypothesis and 
the low density limits. The Michaelis-Menten-Stueckelberg approach substitutes low density of all 
components by low density of the elementary events (or of the correspondent compounds) together with 
the QE assumption. 

The general kinetic law has a simple form: For an elementary reaction 

i i 

the reaction rate is r = cpD,, where > is the Boltzmann factor, f2 = exp (^\ a,/^), fa = —dS/dNi 
is the chemical potential /i divided by RT, and if > is the kinetic factor. Kinetic factors for different 
reactions should satisfy some conditions. Two of them are connected to the basic physics: 

• The detailed balance: The kinetic factors for mutually reverse reaction should coincide, tp + = (p~. 
This identity is proven for systems with microreversibility (Section 7.5). 
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• The complex balance: The sum of the kinetic factors for all elementary reactions of the 
form otiAi —> ... is equal to the sum of the kinetic factors for all elementary reactions 
of the form ... — > J2i a iM (91). This identity is proven for all systems under the 
Michaelis-Menten-Stueckelberg assumptions about existence of intermediate compounds which 
are in fast equilibria with other components and are present in small amounts. 

For the general kinetic law we studied several sufficient conditions of accordance between 
thermodynamics and kinetics: Detailed balance, complex balance and G-inequality. 

In the practice of modeling, a kinetic model may, initially, do not respect thermodynamic conditions. 
For these cases, we solved the problem of whether it is possible to add a linear function to entropy in 
order to provide agreement with the given kinetic model and deformed thermodynamics. The answer 
is constructive (Section 7) and allows us to prove the general algebraic conditions for the detailed and 
complex balance. 

Finally, we have to mention that Michaelis, Menten and Stueckelberg did not prove their "big 
theorem". Michaelis and Menten did not recognize that their beautiful result of mass action law 
produced from the equilibrium relations between substrates and compounds, the assumption about 
smallness of compound concentrations and the natural hypothesis about linearity of compound kinetics 
is a general theorem. Stueckelberg had much more and fully recognized that his approach decouples the 
Boltzmann //-theorem and the microreversibility (detailed balance). This is important because for every 
professional in theoretical physics it is obvious that the microreversibility cannot be important necessary 
condition for the //-theorem. Entropy production should be positive without any relation to detailed 
balance (the proof of the //-theorem for systems with detailed balance is much simpler but it does not 
matter: Just the Markov microkinetics is sufficient for it). Nevertheless, Stueckelberg did not produce 
the generalized mass action law and did not analyze the general kinetic equation. Later, Horn, Jackson 
and Feinberg approached the complex balance conditions again and studied the generalized mass action 
law but had no significant interest in the microscopic assumptions behind these properties. Therefore, 
this paper is the first publication of the Michaelis-Menten-Stueckelberg theorem. 
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Appendix 

1. Quasiequilibrium Approximation 

1.1. Quasiequilibrium Manifold 

Let us consider a system in a domain U of a real vector space E given by differential equations 

^7 = F(aO (112) 
at 

We assume that for any xq E U, solution x(t; xq) to the initial problem x(0) = x$ for (112) exists for all 
t > and belongs to U. Shifts in time, x H- x(t; x ) (t > 0), form a semigroup in U. 

We do not specify the space E here. In general, it may be any Banach or even more general space. 
For nonlinear operators we will use the Freshet differentials: For an operator ty(x) the differential at 
point x is a linear operator (D^l) x : 

d^(x + ay) 



a=0 



dy 

We use also notation (D x ^>) when it is necessary to stress the choice of independent variable. The choice 
of variables is not obvious. 

The QE approximation for (112) uses two basic entities: entropy and slow variables. 

Entropy S is a concave Lyapunov function with non-degenerated Hessian for (3) which increases 
in time: 

^>0 (113) 

dt 

In this approach, the increase of the entropy in time is exploited (the Second Law in the form (4)). 

Formally, any Lyapunov function may be used. Nevertheless, most of famous entropies, like the 
relative Boltzmann-Gibbs-Shannon entropy, the Renyi entropy, the Burg entropy, the Cressie-Read and 
the Tsallis entropies could be defined as universal Lyapunov functions for Markov chains which satisfy 
some natural additivity conditions [32] . 

"Universal" means that they do not depend on kinetic coefficients directly but only on the equilibrium 
point. The "natural additivity conditions" require that these entropies can be represented by sums (or 
integrals) over states maybe after some monotonic transformation of the entropy scale, and, at the same 
time, are additive with respect to the joining of statistically independent systems (maybe, after some 
monotonic rescaling as well). 

Slow variables M are defined as some differentiable functions of variables x: M = m(x). We 
use notation E M for the space of slow variables, M E E M . Selection of the slow variables implies a 
hypothesis about separation of fast and slow motion. In its strongest form it consists of two assumptions: 
The slaving assumption and the assumption of small fast-slow projection. 

The slaving assumption. For any admissible initial state xq E U after some relatively small time 
r {initial layer), solution x(t; Xo) becomes a function of M (up to a given accuracy e) and can be 
represented in a slaving form: 

x(t) = x* M[t) + S(t) for t>r, where M(t) = m(x(t)), \\6(t)\\ < e (114) 
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This means that everything is a function of slow variables, after some initial time and up to a 
given accuracy. 

The smallness of r is essential. If there is no restriction on r then every globally stable system will 
satisfy this assumption because after some time it will arrive into a small vicinity of equilibrium. 

The second assumption requires that the slow variables (almost) do not change during the fast motion: 
During the initial layer r, the state x can change significantly because of fast motion, but the change in 
M = m(x) during r are small with r: The assumption of small fast-slow projection. 

The QE approximation defines the functions x* M as solutions to the following MaxEnt 
optimization problem: 

S(x) — > max subject to m(x) = M (115) 

The reasoning behind this approximation is simple: During the fast initial layer motion, entropy increases 
and M almost does not change. Therefore, it is natural to assume (and even to prove using smallness of r 
and e if the entropy gradient in fast directions is separated from zero) that x* M in (114) is close to solution 
to the MaxEnt optimization problem (115). Further, x* M denotes a solution to the MaxEnt problem. 

Some additional conditions on m and S are needed for the regularity of the dependence x* M on M. It 
is more convenient to discuss these conditions separately for more specific systems. In general settings, 
let us just assume that for given S and m the dependencies m(x) and x* M are differentiable. For their 
differentials we use the notations (Dm) x and {Dx* m )m- The differentials are linear operators: (Dm) x : 
E ->■ E M and {Dx* M ) M ■ E M -> E. 

A solution to (115), x* M , is the QE state, the corresponding value of the entropy 



is the QE entropy and the equation 



represents the QE dynamics. 



S*(M) = S(x* M ) (116) 



™ = (Dm^jFfa)) (117) 



Remark. The strong form of the slaving assumption, "everything becomes a function of the slow 
variables", is too strong for practical needs. In practice, we need just to have a "good" dependence 
on M for the time derivative dM/dt. Moreover, the short-time fluctuations ofdM/dt do not affect the 
dependence M(t) too much, and only the average values 

, wv / \ 1 f t+6 dM 

( M ^ = eJ t -df 

for sufficiently small time scale 9 are important. 



1.2. Preservation of Entropy Production 

Theorem about preservation of entropy production. Let us calculate dS*{M)/dt at point M 
according to the QE dynamics (117) and find dS(x)/dt at point x = x* M due to the initial system (3). 
The results always coincide: 

dS*(M) dS(x) 

dt ~ ~dT (118) 
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The left hand side in (118) is computed due to the QE approximation (117) and the right hand side 
corresponds to the initial system (112). Here, this theorem is formulated in more general setting than 
in Section 2: The slow variables M = m(x) may be nonlinear functions of x. For more details about 
QE approximation with nonlinear dependencies M = m(x) we refer to papers [38,39,66]. The general 
theorem about preservation of entropy production and thermodynamic projector is presented in [39]. 

Proof. To prove this identity let us mention that 

- (DS*) M = (DS*) M o (Hs, (F(x* M )) (119) 



dt v ' V dt / 

where o stands for superposition. On the other hand, just from the definitions of the differential and of 
the time derivative of a function due to a system of differential equations, we get 

dS(x) 



dt 

To finalize the proof, we need an identity 



(DS) x (F(x)) (120) 



(DS*)m o (Dm) x * M = (DS) X * M (121) 
Let us use the Lagrange multipliers representation of the MaxEnt problem: 

(DS) X = A M o (Dm) x , m(x) = M (122) 

This system of two equations has two unknowns: The vector of state x and the linear functional A M on 
the space of slow variables (the Lagrange multiplier), which depends on M as on a parameter. 
By differentiation of the second equation m{x) = M, we get an identity 

(Dm) x * M o (Dx* M ) M = id EM (123) 

where id is the unit operator. 

Lagrange multiplier Am is the differential of the QE entropy: 

(DS*) M = A M (124) 

Indeed, due to the chain rule, (DS*) M = (DS) X * M o (Dx* M ) M , due to (122), (DS) X = A M (Dm) x 
and, finally 

{DS*) M ={DS) x * m o (Dx* m ) m = A M o (Dm) x o (Dx* M ) M 
=A M o id EM = A M 
Now we can prove the identity (121): 

(DS*) M o (Dm) x * M = A M ° (Dm) xh = (DS) X * M 

(here we use the Lagrange multiplier form (122) again). □ 

The preservation of the entropy production leads to the preservation of the type of dynamics: If for the 
initial system (112) entropy production is non-negative, dS/dt > 0, then for the QE approximation (1 17) 
the production of the QE entropy is also non-negative, dS*/dt > 0. 

In addition, if for the initial system (dS/dt) x = if and only if F(x) = then the same property 
holds in the QE approximation. 
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2. First Order Kinetics and Markov Chains 

First-order kinetics form the simplest and well-studied class of kinetic systems. It includes 
the continuous-time Markov chains (the master equation [67]), kinetics of monomolecular and 
pseudomonomolecular reactions [17], and has many other applications. 

We consider a general network of linear reactions. This network is represented as a directed graph 
(digraph) ([33,68]): Vertices correspond to components Bj (1 < j < q), edges correspond to reactions 
Bj — > B[. For each vertex Bj a positive real variable qj (concentration) is defined. For each reaction 
Bj — > B\ a rate constant «y > is given. To follow the standard notation of the matrix multiplication, 
the order of indexes in Kji is always inverse with respect to reaction: It is Kj^, where the arrow shows 
the direction of the reaction. The kinetic equations for concentrations 9 have the form 

7j| = Y ~ K &j) ( 125 ) 
The linear conservation law (for the Markov chains this is the conservation of the total probability) is: 

^ = i.e., J2 - K ifr) = • (126) 

Let a positive vector q* (q* > 0) be an equilibrium for the system (125): For every j = 1, . . . , q 

Y = Y K &j ( 12V ) 

l.lft 

An equivalent form of (125) is convenient. Let us use the equilibrium condition (127) and write 



Y = Y «W )% = Y K fl^% 



Therefore, under condition (127) the master equation (125) has the equivalent form: 

The following theorem [69] describes the large class of the Lyapunov functions for the first order 
kinetics. Let h(x) be a smooth convex function on the positive real axis. A Csiszar-Morimoto function 
H h (q) is (see the review [32]): 



The Morimoto //-theorem. The time derivative of Hh(q) due to (125) under condition (127) 
is nonpositive: 



dH h k) 
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Proof. Let us mention that for any q numbers h { , j^i K ij^j(hj ~ hi) = 0. Indeed, for hi = Pi/p* 
this is precisely the condition of conservation of the total probability for equations (128). The extension 
from a simplex of the hi values (J2iPi^i = L hi > 0) to the positive orthant M. 9 + is trivial because of 
uniformity of the identity. Finally, if a linear identity holds in a positive orthant then it holds in the whole 
space M. q . Therefore, 



< 



The last inequality holds because of the convexity of h(x): h'(x)(y — x) < h(y) — h(x) 
(Jensen's inequality). □ 

For example, for the convex function h(x) = x(lnx — 1) the Csiszar-Morimoto function is: 

^w = E^( ln (|) - 1 ) ( 13 °) 

This expression coincides with the perfect component of the free energy (50) (to be more precise, 
/ = RTH h (q)). 

Each positive equilibrium <;* belongs to the ray of positive equilibria A<^* (A > 0). We can select a 
/x-normalized direction vector and write for an equilibrium ^ cq from this ray: 

? 6q = 13; 



where /3 = ^. S eq . 

The kinetic equations (125) allow one and only one ray of positive equilibria if and only if the digraph 
of reactions is strongly connected: It is possible to reach any vertex starting from any other vertex by 
traversing edges in the directions in which they point. Such continuous-time Markov chains are called 
ergodic chains [70]. 

Let us assume that the system is weakly reversible: For any two vertices Bi and Bj, the existence 
of an oriented path from Bi to Bj implies the existence of an oriented path from Bj to Under this 
assumption the graph of reactions is a unit of strongly connected subgraphs without connections between 
them. Let the graph of reactions Bj — > Bi have d strongly connected components C s and let V s be the 
set of indexes of those Bj which belong to C s : Bj E C s if and only if j E V s . For each s = 1, . . . , d 
there exists a conservation law 

0s(s) = ^ 9 = cori ^ (131) 

For any set of positive values of /3 S > (s — 1, . . . , d) there exists a unique equilibrium of (125), <^ eq , 
which is positive (<^ cq > 0). This equilibrium can be expressed through any positive equilibrium ^*: 

c cq = R ^lll (132) 

For positive (3 S , logarithms of <^ eq form a c/-dimensional linear manifold in R q : 
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The natural coordinates on this manifold are In f3 s . 

Remark. In the construction of the free energy (130) any positive equilibrium state ? eq can be used. 
The correspondent functions differs in an additive constant. Let us calculate the difference between two 
"free energies" : 




(134) 



The result is constant in time for the solutions of the master equation (125), hence, these functions are 
equivalent: They both are the Lyapunov functions for (125) and have the same conditional minimizers 
for given values of f3 s > (s — 1, . . . , d). 
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